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ABSTRACT 


The three dimensional transient temperature variations during autogenous Gas 
Tungsten Arc Welding are determined. The model emplovs a combination of un-equally 
spaced moving meshes to minimize the total number of nodes. Finite differencing is used 
for the spatial terms. The resulting ordinary differential equations for the transient ev- 
olution of thermal transport are solved using the fourth order Runge-Kutta technique. 
The temperature dependent thermal properties and latent heats of phase transformations 
are accounted for. Computations are carried out for a rectangular parallelepiped, with 
convective and radiative surface thermal conditions. Results are presented for the evo- 
lution of thermal profiles during ideal welding conditions. These are next compared with 
variations obtained due to defects such as weld track misalignment and inclusions. A 
study of the startup and shutdown transients makes it possible to control the cooling 
rate during these transients. The potential use of this model in the development of an 
expert Welding system using infrared imagery is indicated. In addition, a low cost infra- 
red detector using an indium-arsenide diode is prototvped to determine its feasibility for 


production welding control. 
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I. INTRODUCTION 


The Navy's David Taylor Research Laboratory is currently in the process of devel- 
oping an expert welding svstem. This integrated automatic welding svstem will be used 
for the Gas Metal Arc (GMA) welding of submarine hulls. This technology is mandated 
by the problems in mass production welding with the new high strength, low allov steels. 
These steels have very tight limits on cooling rates which result in very low metal depo- 
sition rates. To meet production requirements, the use of an expert welding system is 
deemed mandatory. One portion of the expert welder is the infrared vision system. To 
effectively use this svstem it is desired to determine how various welding parameters af- 
fect the thermal signature at the surface of the plate being welded. However, no numeric 
models have been developed to supplement these experiments. Due to the cost and time 
associated with laboratory experiments, a numeric model of the welding process was 
needed to aid in this development. 

Three specific problems were identified as candidates for investigation of their effect 
on the surface temperature profile. These were: 

e The effect of inclusions in the weld. 


e The effect of weld seam-to-arc misalignment 


e The effect of start-up and shut-down on cooling rates. 


The first two of these, by definition, are unsymmetrical problems. This precludes the 
use of symmetry in the computer model, something most models use to control the 
number of nodes to solve. In general, the solution cost is at best based on the square 
of the number of nodes for anv solution technique requiring matrix inversion. However, 
since accurate results are desired, a three dimensional, variable coefficient model 1s the 
requirement. 

The use of the finite element method for three dimensional weld modeling has been 
under development for a number of vears. The current finite element models have in- 
cluded features such as temperature dependant material properties, complex and chang- 
ing geometries and transient heat conduction analysis. The foremost of these is Goldak 
{Ref. 1] who concludes that present day computers are not capable of solving meaningful 
three dimensional problems. He estimates that it would take a CRAY-2 computer about 


6 hours to solve a 10 second problem of interest using the finite element method with 


only 659 nodes. This is because the sparse matrix resulting from the formulation must 
be solved using an iterative process such as Gauss-Seidel. Two researchers, Mangonon 
and Mahimkar [Ref. 2], have successfully developed a simple three dimensional model. 
When using 2500 nodes they took four hours of CPU time to simulate six seconds of real 
time. They used a grid spacing of 2.12 millimeters, which allowed them to take larger 
time steps, instead of the 1.0 millimeter grid spacing recommended by Goldak. To re- 
duce the number of nodes, most researchers have invoked svmmetry whenever it was 
possible. The end result of these difficulties is that welding models have been limited in 
scope. Most are limited to predictions of the weld pool geometry and the cooling rate 
in the heat affected zone, since these are easv to obtain. 

These computational limitations, compounded by the need of a model that does not 
use symmetry, forced a reconsideration of the modeling technique to use. The explicit 
finite difference technique is the simplest and most direct of all of the numeric techniques 
for solving the heat conduction equation. Since the equations are not coupled, no ma- 
triX inversion 1s required to solve the system, greatly reducing the computational effort 
associated with large numbers of nodes. However, there remains the classic limitation 
Onethe line steprone om +. Where Fo is the Fourier number, and the 6 1s for a the. 
dimensional model. Of course, the Fourier number is not the only limitation on time 
step. The degree of accuracy is also dependant on the time step. Because of the large 
temperature gradients in the vicinity of the arc, very small time steps are required for 
accuracy. Unlike most time dependant problems, this condition is not transient, and 
dynamically changing the time step is not an option. Testing with a simple two dimen- 
sional model showed that for the temperature gradients of interest, the time step required 
to be within one percent of the time step independent solution was more limiting than 
the Fourier number limitation. Further numeric testing showed that the equivalent im- 
plicit finite difference model was just as accurate (or inaccurate) as the explicit finite 
difference model under these conditions, while having the more expensive cost of matrix 
solving. 

With this insight into the unique problems of welding analysis, that accuracy con- 
trolled the time step, not the stability limit of the Fourier number, the explicit technique 
becomes a viable alternative. To allow taking as large a time step as possible and still 
have reasonable accuracy, the finite difference equations are written in the semi-discrete 
form. This allowed applving fourth order Runge-Kutta. The combination of the explicit 


technique and Runge-Kutta proved to be the key to the solution. This resulted in a very 


to 


fast, accurate three dimensional model of the arc welding process which was able to 
handle grids with large numbers of nodes. 

Using these techniques the model was wntten and debugged. A series of test were 
conducted to compare the results of the model with the available analytical solutions and 
experimental data. Following the verification, each of the major goals was investigated. 
In addition, some simple laboratory experiments were designed to verifv the specific re- 
sults of the model. The model proved to be very versatile in answering the posed 
questions and is a valuable tool for continued research. 

In the course of this research, a low cost alternative to the infrared camera system 
was proposed. In lieu of the expensive and delicate vision system, the use of a simple 
infrared diode was considered. There exists a commercially available indium arsenide 
diode which has an acceptable detection range. Experimental work was performed to 
determine if the diode had the sensitivity and accuracy necessary to be a useful industrial 
detector. The results of this testing, especially in light of the results of studving flaws 
and arc-to-seam misalignments, indicates that a simplified vision system would be a cost 


effective option. 


Il. THE THREE DIMENSIONAL WELD MODEL 


A. SELECTION OF MODELING TECHNIQUE 

Both the implicit and explicit finite difference techniques were considered. The 
benefits of the explicit technique are many. The code is easy to write for the three di- 
mensional case since node ordering 1s not important. The memory storage requirements 
are substantially less since an N by N matrix 1s not being created. The solution cost is 
alwavs N calculations per time step. The formulation can be developed in the semi- 
discrete form, where only the spatial derivatives are written using difference equations. 
The time derivative is then evaluated using fourth order Runge-Kutta. For a rectangular 
Steel plate, the boundary conditions are simple to encode. Another possible technique, 
the implicit method is not limited by stability in choosing a maximum allowable time 
step based on the Fourier number. This time step limitation is discussed in detail in [Ref 
3]. The benefit of the finite element method is that it easily adapted to the weld pool ge- 
ometrv and it 1s simpler to increase the nodal density in the vicinity of the arc. 

Because of the simplicity of the explicit finite difference technique, a series of tests 
were performed to compare the iumnplicit technique to the explicit technique formar 
problem of interest. These tests compared both methods, solving a simple problem with 
an analytical solution. The methods were then applied to a problem of a step transient. 
It was found that with the large temperature gradient of a step, the limitation on the 
Fourier number was not what limited the time step size. For a fixed grid size, it was 
found that both the implicit and explicit methods developed a greater than three percent 
error When the time step was increased to the Fourier number critical time step. The 
error 1s measured relative to a time step independent solution, which is obtained by 
halving the time step until the solution changes by less than 0.1 % . This solution 1s 
used in lieu of an analytical solution to calculate the error. In addition, the explicit 
method was always the more accurate of the two methods when solving problems with 
large temperature gradients. Thus no benefit would be gained bv using the implicit 
technique. 

To keep the number of nodes down to a solvable number but still have a fine grid 
mesh in the region of large temperature gradients, a combination of three mesh sizes is 
used. These meshes shown in Figure | include the fine mesh, medium mesh and coarse 


mesh. Each small rectangular region corresponds to a control volume with the 
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The Coarse, Medium and Fine Meshes 


Figure 1. 


associated node located at the centroid. For the case of surface nodes, the node 1s ac- 
tuallv at the surface and not at the centroid of the contro! volume; for the corner nodes 
of surfaces, the node is at the corner. This is consistent with the normal method of as- 
signing nodes for finite differences. In the figure, this 1s shown by the difference in 
control volume size, since for a given mesh, the nodal spacing is constant. 

The spacing of the fine mesh is one millimeter, the medium mesh 1s three millimeters 
and the coarse mesh 1s nine millimeters. The factor of three 1s chosen to allow smoother 
shifting between different sized meshes. The fine mesh moves with the arc through the 
medium mesh in three millimeter jumps. After the arc moves nine millimeters, the me- 
dium mesh shifts nine millimeters in the coarse mesh and the fine mesh 1s back at its 
initial point with respect to the medium mesh. In this wav the arc is always located in 
the fine mesh grid while the size of the fine mesh grid is Kept to a minimum. 

To allow taking as large a time step as possible and still have good accuracy, the 
finite difference equations are written in the semi-discrete form. This allowed applying 
fourth order Runge-Kutta. Not only does this allow taking time steps up to the Fourier 
number limit, but it also stabilizes the system of equations and it is possible to take time 
steps in excess of the Fourier number limit and still have reasonable accuracy. Fora 
Fourier number limitation of eight milliseconds, a time step of ten milliseconds is rou- 
tinelv used with an error of less that one percent. This error is based on taking smaller 
time steps until the solution converged. As before, this converged solution 1s the refer- 


ence for the error. 


B. CREATING THE MESHES 

The geometry selected 1s a thick plate of steel. The selected size to model/1Swiaaen 
7” by 12” steel. The model actually uses 27mm by 180mm bv 315mm. The whole plate 
is modeled by the coarse zone using a two dimensional grid 21 by 36. The medium grid 
is centered about the arc and is 27mm by 81mm by 81mm for a grid of 10 by 27 by 27. 
The fine grid is centered about the arc, but does not penetrate the plate, and covers a 
size of 7.5mm by 27mm by 27mm for a grid of 8 by 27 bv 27. This gave a total of 13,878 
nodes. The coarse zone 1s made two dimensional since it is far enough away from the 
arc for the temperature profile to be essentially isothermal in the z (vertical) direction. 

Pointers are used to keep track of the relative positions of the grids. When the grids 
are moved relative to one another the temperature of the finer grid nodes being left be- 
hind are averaged to provide a temperature to the coarser node which would occupy the 


same space. The finer grid nodes moving to where the coarser grid had been are all 


assigned the temperature of the coarser grid point whose space they now occupied. The 
nodes are shifted only in increments of three of the finer grid, which corresponds to one 
increment of the coarser grid. 

The problem of temperature dependant material properties, including phase trans- 
formation, 1s also treated. The material properties are constant in both the coarse and 
medium mesh zones and are evaluated at room temperature. In the fine zone, the 
change in enthalpy is calculated instead of the change in temperature. The correspond- 
ing temperature for a given enthalpy is found using a piece-wise continuous function. 
This took into account the austenite to ferrite transitions as well as melting. The tem- 
perature dependence of thermal conductivity is also similarly handled. The sizing of the 
fine zone insured that the high temperatures where these effects dominated are inside the 
fine zone mesh. In addition, a radiation boundary condition is added in the fine zone. 

The use of different sized grids creates an additional problem of calculating the heat 
transfer at the interface zones. The bookkeeping becomes involved and a detailed ex- 
planation of how the program is written and the source listing are included in appendix 


A. A brief discussion of how the equations are set up on the mesh follows. 


C. FINITE DIFFERENCE DERIVATION 
1. The Fine Mesh Zone 

The fine mesh region uses the enthalpy form of the finite difference conduction 
heat equation. In addition, the thermal conductivity 1s a function of temperature and 
so will be different at each node. The thermal conductivity used for each difference 1s 
the harmonic mean for the two nodes being differenced. A side view of the fine mesh and 
its interface with the medium mesh 1s seen in Figure 2. Some of the representative en- 
ergy balance equations are listed next. 

a. Interior Nodes 


The most general case is an interior node, given in equation 2.1. 
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Side View of Fine and Medium Interface 


Figure 2. 


For programming, equation 2.1 can be simplified by the function defined in equation 2.2. 


This calculates the harmonic mean of the thermal conductivity at each nodal pair. 
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Soo) — =, ; (2.2) 
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b. Surface Nodes 
For surface nodes there are three terms that contribute to the change in 


enthalpy. As shown in equation 2.4 these are the conduction, convection and radiation. 


AH = ON ee, a5 Gaon a Grad (2.4) 


The convection and radiation terms are given by. 





} 
M4eony,., = 7 ae Tiyx) (2) 
ee eee 2.6 
Mrad,,, = Az ( inf ales (2.6) 


For the surface nodes, the node immediately below the surface has twice the effect of the 
nodes on the side, because it has twice the heat transfer area. Thus @,,,, must be modi- 
fied accordingly. 

2. The Medium Zone Mesh. 

In the medium zone, all of the material properties are held constant. In this 
case it 1s simpler to write the equations in the temperature form of the finite difference 
equations. The material properties are evaluated at 300 k. The relationship of the 
medium mesh to the coarse mesh 1s shown in Figure 3. 

a. Interior Nodes 
Equation 2.7 shows the tvpical form for an interior node, where Fo 1s the 
Fourier number. 
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Figure 3. Side View of Medium and Coarse Interface 


The medium zone is one of the more complicated to program since it interfaces with 
both the fine and the coarse grids. Interfacing problems are discussed in “4. The 
Interfaces” on page Il For the surface nodes an additional term 1s added to take into 
account the free convection. Care is taken to note that the surface nodes have only 
one-half of the volume of an interior node. Equation 2.8 shows a typical form for a 


surface node, 


n+l - ae ip : 
AT, 5; = FOU say ia Dail eee tee eee lee, 
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where Biis the Biot number for free convection. 
3. The Coarse Zone Mesh 
In the coarse zone, all of the material properties are constant. Since this models 
the plate far from the arc, it is only two dimensional in the X-Y plane. These equations 
are also written in the temperature form of the finite difference equations. The top and 
bottom surface are free convection while the ends are considered adiabatic, since the 
plate size for the model is based on selecting a large enough area for this to be true. 


Equation 2.9 shows the typical form for a coarse mesh node. 


ATS = FolTiyy + Tey + Tiger + Tye + 2BiT iy — (4 + 2Bi)T},] (2.9) 


Care is taken with the top and bottom Biot numbers since the heat transfer area at a 
node's surface is one third that of the node’s sides. 
4. The Interfaces 

The interfaces are handled by treating them as boundaries. Each grid mesh is 
solved separately, using the temperature from the boundaries of adjacent meshes. 
Pointers are used to keep track of the relative position of the fine mesh to the medium 
mesh and the medium mesh to the coarse mesh, since these move. Because the control 
volumes for the nodes are different at the interface, the finite difference equation is mo- 
dified. A one dimensional example is shown to demonstrate this. For the case where 
nodal point m-1 is twice as far from nodal point m as nodal point m+ 1, the second de- 


rivative may be expressed by: 
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These first derivatives may be in turn expressed as: 
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These combine to form the following approximation to the second derivative at nodal 
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a. The Fine Zone to Medium Zone Interface 
The details of the medium to fine interface are shown in Figure 2 and Fig- 
ure 4. The fine to medium zone interfaces have two complicating features. The first 1s 
that the fine zone uses variable thermal properties while the medium zone uses constant. 
To insure conservation of energy, the constant thermal properties of the medium zone 
are used when calculating the change in enthalpy at the interface. This is done because 
using the same thermal properties insures that the amount of heat entering the medium 


zone is equal to the amount of heat leaving the fine zone. The second complication 1s 
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The Interface Between the Medium and Fine Meshes 


Figure 4. 
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This is seen in the last term of Equation 2.10 


the difference in node control volumes. 


‘-k’ direction. 
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the interface relation for a non-surface, 
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The factor of 1 21s due to the medium node being twice as far away as a fine node, with 
the change in enthalpy being related linearly to the diffusion distance. The heat transfer 
area for a medium node to one fine node is the same as a fine node to one fine node, so 
it does not effect the result. Figure 2 shows these relationships. 
b. The Medium Zone to the Fine Zone 

Interfacing the medium grid to the fine grid is a two step process. First, for 
everv node of the medium interface, there are nine fine nodes at the interface. This is 
handled by taking the average of the nine nodal temperatures. This average temperature 
is then considered to be at a pseudo nodal point, which is only 2'3 the distance away as 


a normal medium nodal point. The energy balance equation then becomes, 
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As can be seen in equation 2.11, this adds a factor of 3,2 to the equation. A quick cal- 
culation will show that energy is conserved; the heat flux calculated by the fine mesh is 
equal to the heat flux calculated by the medium mesh. 
c. The Medium Zone to the Coarse Zone 
The interface of the medium mesh to the coarse mesh has identical geometry 


as going from the fine mesh to the medium mesh. The energy balance equation is, 
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As can be seen in equation 2.12, the factor of 1/2 is present for the same reason as in the 
fine to medium interface. It should be pointed out that the Fourier number for each 
mesh is different due to the difference in nodal spacing. The Fourier number in an 
equation is always calculated for the spacing of the nodal point being evaluated. 
d. The Coarse Zone to the Afedium Zone 
The coarse to medium interface is handled in the same manner as the me- 
dium to fine interface. There are 30 adjacent medium nodes to one coarse interface node 


and these nodal temperature are averaged together: 
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Equation 2.13 1s the energy balance equation for this case. Again, the pseudo node is 
only 2.3 the distance awav as a normal coarse node and the factor 3/2 appears. 


Equation 2.13 shows the relationship at the interface in the ’-}’ direction. 


D. SIMULATING THE ARC 

Simulating the arc and the weld pool is one of the most difficult aspects of the mo- 
del. Manv researchers shape the arc to obtain the desired empirical results of weld pool 
size, shape and penetration. Such techniques, however, are not appropriate for the study 
of transients for they only hold in the steady-state. It is desired to add heat to the weld 
with little or no shaping and yet still model known weld pool dynamics. Two different 
approaches are taken, both having appropriate applications. 

In the first approach the arc and the weld pool are simulated by spatially shaping 
the forcing function. Based on the experiments of Lu [Ref. 4] the energy of the arc can 
be assumed Gaussian in distribution. This is then added volumetrically to the fine mesh. 
The arc is henuspherical with a radius of four millimeters. Further shaping of the arc 
by use of double ellipsoids as discussed by Goldak [Ref. 5] are not used because it in- 
terfered with the study of transients. Since this weld pool simulation shapes the weld 
pool it is only useful for studying transients that do not pass through the weld pool nor 
effect weld penetration. The weld pools generated are shallow, with little penetration. 
This 1s tvpical of reverse circulation weld pools where flow is from the center of the pool 
surface to the edge of the pool. 

In the second approach, the arc 1s given a Gaussian energy distribution in the sur- 
face plane but the energy 1s added only to surface nodes. To simulate the penetrating 
effect of the arc pressure, the thermal conductivity is made directional. Vertical thermal 
conductivity 1s an order of magnitude greater than the horizontal. The directional ther- 
mal conductivity only applies in the molten metal temperature range. The weld pools 
generated showed deeper penetration, which is typical of normal circulation, where the 
flow is from the center to the bottom of the weld pool. 

In both cases the program allowed for the arc to move in the X-Y surface plane. 
The amount of heat added to each node 1s determined by calculating the Gaussian dis- 
tribution factor for each node and summing them. This provided a normalizing volume 
by which to divide the arc input energy for that time increment. The normalized arc 


power density is then used to calculate the enthalpy change of each node. 
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FE. APPLYING RUNGE-KUTTA 

The application of Runge-Kutta to the explicit finite difference problem proved to 
be a very powerful tool. Many difficult problems have been successfully solved using 
explicit Runge-Kutta as described in [Ref. 6]. In appendix B it is shown that for the 
constant coefficient, three dimensional case the limitation imposed by the Fourier num- 
ber is eliminated altogether. Thus it is possible to take larger time steps and still main- 
tain accuracy. However, due to the large temperature gradients in the vicinity of the 
heat source, it is found that time steps in excess of the Fourier number limitation had a 
rapid degradation in the accuracy and stability of the solution. 

The general solution technique of Runge-Kutta involves solving the equations four 
times for each time step. The advantage is that this provides an accuracy of O(h‘), a 
substantial improvement over the O(h') of using forward differencing. Numeric testing 
showed that about 1’20 the time step would be required to obtain the same accuracv 
(even with implicit), so that Runge-Kutta resulted in a times savings of at least 80 per- 


Celt. 


F. PROGRAM OVERVIEW 

The program sequence 1s as follows. The first segment either defines the initial 
conditions or reloads a previous problem. The main processing portion consists of three 
parts. First the heat of the arc is added to the fine mesh. Second the change in tem- 
perature is calculated by applying Runge-Kutta. Third, the mesh is checked to see if it 
is time to move meshes. If it is, the meshes are moved and the pointer updated. The last 
two parts are associated with output. One section handles running output during exe- 
cution and the last portion saves the final temperature distributions and data necessarv 
to restart the problem, if desired. 

There are three major subroutines in the program. These are labeled FIN. MED 
and COR for fine, medium and coarse and are where the finite difference equations are 
set up. These are called by the Runge-Kutta solver. There are also three functions. 
One finds the thermal conductivity as a function of temperature and the other two are 


used for converting to and from enthalpy and temperature. 


IH. VALIDATION OF THE NUMERIC MODEL 


The accuracy of the model in reflecting the actual welding conditions is heavily de- 
pendant on the assumptions made in stating the problem, the accuracy of the material 
properties and the resolution of the numeric technique. In addition, it is desirable to 
Validate the results of the model against other models and analytical solutions were they 
are available and against experimental results. This chapter will discuss each of these 


areas in some detail. 


A. VALIDATION OF THE MODEL ASSUMPTIONS 

The assumptions of the model are discussed in detail in each of the following sec- 
tions. The decisions leading to these assumptions were based on the work of previous 
investigators, Goldak {Ref. 5: pp. 587-600] having provided the best summarv in his ar- 
ticle. In addition, several simple two dimensional models were used to confirm the ac- 
curacy and stability of several solution techniques to further refine the approximation. 
There are two conflicting desires in developing any model, the first is to maximize the 
resolution of the solution and the second is to minimize the computational effort. To 
achieve this, it is necessary to insure that greater resolution is applied only were it is 
Warranted. 

Perhaps the first validation of the model is to look at the actual results of the model. 
Figure 5 shows the surface temperature variation for the x-v plane of the fine zone 
during a tvpical startup, Figure 6 is a contour plot of the same startup, and Figure 7 is 
a contour plot of the x-z plane directly under the arc. These all show the expected 
shapes of the weld pool and thermal contours during a startup sequence. In the surface 
profiles, the portion of the Weld pool not under the arc 1s seen as a region of near con- 
Stant temperature due to the phase change occurring. In all of these pictures the arc 1s 
at av position of about 17 mm. 

1. Fine mesh spacing of one millimeter provides adequate resolution 

This is the recommended spacing of Goldak, though other researchers have 
successfully used larger spacings. Referring back to Figure 5, it can be seen that this 
spacing in the fine zone is more than adequate to produce a smooth thermal contour 
of the temperature. However, in those cases were it is desired to measure a temperature 
difference from the quasi-steadystate value, such as flaw detection, some model noise 


Was notice that could have been limited bv decreasing the mesh spacing. 
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Surface Temperature Variation in the x-y Plane: The torch moves 
parallel to the fine y axis atx = 14mm. Q = 2544 watts, v = 4 mm/s. 


Figure 5. 
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Figure 6. Surface Temperature Contours in the x-y Plane: The torch moves 
parallel to the y axis at x = 14 mm. The x-y coordinate system is fixed 
to the fine zone. Q = 2544 watts, v = 4 mm/s. 


2. Coarse and Medium meshes decrease the number of nodes 
The idea of using different mesh sizes is to minimize the the number of 
nodes[Ref. 7]. Unlike the Finite Element Method (FEM) it is not a good practice to 


gradually change the mesh size for Finite Difference. This is because the accuracy of the 
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Finite Difference Method is O(/17?) for uniform mesh sizes versus O(/7) for variable mesh 
sizes. It is necessary to insure that the gradients are low when shifting to a coarser mesh 
tO minimize the error induced by the volumetric averaging. This is checked visually by 
studving the temperature profiles that the model generates. For instance, in the case of 
studying cooling rates, it was found that the region behind the arc needed to be increased 
to allow the cooling rates to be determined more accurately. This particular modifica- 
tion is discussed in further detail in this chapter and caused the program to take ap- 
proximately twice as long to run. 
3. Radiation could be approximated by a linear coefficient 

This was a simplification used to improve the speed of the program. It was 
based on the fact that for a thick plate, conduction dominates the problem compared to 
radiation and that in this regions, the surface temperature 1s below the temperature at 
which significant heat will be lost in radiation. Any surface temperature losses were si- 
mulated by use of h, the convection term. This was applied only in the medium and 
coarse regions. 

4. The Free Convection Coefficient may be approximated as a constant 

Again, this is a simplification for two reasons. The first 1s that 1t speeds up the 
operation of the program and the second 1s that the value of h in the welding environ- 
ment is not accurately known anyway. In addition, this effect 1s small compared to 
conduction for the thick steel plate being modeled. Work done by Goldak [Ref. 5: pp. 
587-600] confirms the practicality of this assumption. 

5. The energy of the arc may be added as a volumetric gaussian source. 

This is discussed in greater detail in Chapter 2 where it was mentioned that this 
has been demonstrated both by other models and experimental work on are power dis- 
tributions. In this model, no attempt was made to shape the arc power distribution to 
obtain a desired matching with an experimental weld pool. Arc shaping was avoided to 
since the response of the weld pool was one of the desired results and this would be ob- 
scured if shaping was used. 

6. Material properties are constant far from the arc 

This is used in the medium and coarse regions since all of the elevated temper- 
atures are contained within the fine region. Thus, over the lower temperatures that are 
experienced in the these regions, this is a fair approximation. This has a significant im- 
pact On program speed, since it 1s much faster to calculate the finite difference equation 


with constant coefficients. 
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B. INITIAL PROGRAM DEBUGGING 

In writing a program from scratch a large number of programming errors occurred, 
most of which were clerical in nature. These errors were found by letting the program 
run and observing the results, in most cases the mistakes clearly stood out as violations 
of the the laws of thermodymanics. Next it was desired to insure that the arc energy 
was correctly inputted. A simple energy balance was performed to insure that the heat 
from the arc is conserved. After the model passed these simple tests it was considered 


free of basic programming errors. 


C. COMPARISON TO AN ANALYTICAL SOLUTION 

The analytical solution available is called Rosenthal’s Solution, and it solves the 
temperature distribution due to a point source of energy moving along the surface of a 
semi-infinite solid at uniform velocity. Rosenthal [Ref. 8] describes the derivation of this 


solution in detail and the result is presented as equation (3.1). 


( as 
i 1G: = = D (r+x)/2x (3.1) 


where 


ee ae) een” 


and rv is the velocity of the arc. The coordinate svstem is centered at and moves with the 
arc. This elegant solution is valid for only constant thermal properties. Thus, for the 
purposes of validation, the model 1s modified accordingly. All of the thermal properties 
are made constant (thermal conductivity and the heat capacity). The arc power 1s re- 
duced to avoid melting and heat 1s added to a single node instead of using a volumetric 
source. Also, since Rosenthal’s solution does not take into account radiation and con- 
vection, these are eliminated from the model. 

The results of the validation are shown in Figure $8, which is a comparison of 
Rosenthal’s solution to the model. The current arc position 1s y=17, and the temper- 
ature contour plotted is at x=z=0. This demonstrates that the finite difference model 
is converging to the analvtical solution. The difference seen is the expected one, since 
any finite difference model will be stiffer than an analytical one. The largest source of 
error is at the discontinuity at the center of the arc, which 1s expected due to the point 
source. If the model is attempting to simulate a point source, it 1s clear that a finer or 


different mesh scheme would be required near the arc. However, since it is desired to 
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Figure 8. Rosenthal’s Solution Compared to the Finite Difference Model 


model a volumetric source rather than a point source, the present mesh size is adequate. 
This 1s supported by Goldak [Ref. 5: p. 590] who reports that having at least four ele- 
ments in the arc is adequate and eight is preferable. The spacing of | millimeter provides 


a nunimum of eight nodal points under the arc. 


D. QUASI-STEADYSTATE THERMAL PROFILES 

There are two comparisons that can be easily made to experimental data. This is 
the measured cooling rates and the weld pool dimensions. Both of the theoretical and 
experimental work is summarized by Lancaster in his text on welding [Ref. 9]. In gen- 
eral, the analytical model provides results within a factor of 2 of the experimental results. 


To characterize the weld, a non-dimensional parameter, n, is defined by equation (3.2). 


qv 
y= 


oe (3.2) 
dna C,(8, — 9p) 


Average values of the weld parameters are: 


g = 2544 warts 
vy = 0.004 pn/s 
% = 9.0x107° m/s 
C= 4x10° watts|m>K 
0. = reference temperature 


0, = 300 degrees K, initial temperature 


Using the melting point of steel, 1723 degrees Kelvin, as the reference temperature, it is 
found that n is 1.8. Referring to figure 3.17 of Lancaster, the weld pool width is eight 
millimeters, which 1s exactly what 1s predicted by the model. 

As discussed by Lancaster, there is a great deal of data scatter in trying to determine 
cooling rates, and in general, the experimental cooling rates are less than the theoretical. 
This is only in part due to the finite size and time delays associated with thermocouples. 
To study this, a modified version of the welding program called WELDC was used. As 
discussed in Appendix C, this program uses extended fine and medium regions to 1m- 
prove the accuracy of cooling rate measurements, since the temperature at which we 
Wish to measure the cooling rate typically occurs 30 mullimeters behind the arc. 
Figure 9shows the cooling rate at the arc centerline for an arc startup and shut down 
transient. The reference temperature is 810 degrees Kelvin and the initial temperature 


is 300 degrees Kelvin. In this case, the arc power was only 1950 watts so that the 


ZS 


LEGEND 
da DEG C 
AC Oi 


: 
O 
Y) 
ee 
fx] 
Pe 
O 
= 
eas 


Q=1950 WATTS V=4 MM/S 


ce a0) 30 40 00 60 70 
ARC TRAVEL (MM) 





Figure 9. The Cooling Rate Measured at 535 deg C: Shown for an arc startup 


and shutdown sequence, taken along the arc centerline. 


characteristic parameter n was 3.8. Referring to figure 3.19 of Lancaster, the expected 
value is 140 degrees C per second. This is in excellent agreement with the quasi- 
steadystate value shown in Figure 9. The difference of ten percent is easily explained 
due to the error in averaging the material thermal properties in calculating the charac- 


teristic parameter n as well as the experimental data scatter. 


FE. EFFECT OF ARC OSCILLATION 

Most automatic welders use an oscillating arc to improve the weld quality. This is 
because the transverse oscillation insures that the weld pool is spanning both sides of the 
weld seam, preventing a lack of fusion. A trial run was performed to compare the ther- 
mal signature of a steady weld to the case of an oscillating one. For this test the arc 
Was given a sinusoidal oscillation at four hertz and a peak to peak transverse amplitude 
of one millimeter. The results of this are shown in Figure 10 and as expected the oscil- 


lating arc has a wider thermal profile and slightly lower peak temperatures. This again 
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Figure 10. Comparison of a Steady Versus an Oscillating Arc.: a. Steady arc 


temperature contour. b. Oscillating arc temperature contour. 


shows the usefulness of this model, in that it can easily simulate complex arc patterns 


since it does not invoke symmetry in the geometry. 


F. DETERMINING THE MATERIAL PROPERTIES 

All of the initial results of the model are based on using the material properties of 
a low carbon steel provided by Goldak [Ref. 5: p. 591]. However, since the actual test 
pieces were made of HY-80 steel it was desired to use the properties of this steel also. 
The data available on the HY series of steels was provided by Morris (Ref. 10]. This 
data provided thermal conductivity up to 811 degrees K and specific heat up to 673 de- 
grees K. The density of the material 1s considered to be constant and 1s about 
7890kg/ne? The thermal coefficient of linear expansion up to 1100 degrees F is 7.1x10-° 
in./in./F. which is less than a one percent change. Above this there is the phase change 
to face centered cubic where the metal actually contracts, hence considering the density 


constant is a good approximation [Ref. I]: p. 13]. 
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The thermal conductivity of metals in genera] does not have a theoretical basis at 
low temperatures but experimental work has shown that above 1200 degrees K for iron 
based alloys, all alloys have the same thermal conductivity. In Figure 1] is the plot of 
thermal conductivity for pure iron, a stainless steel and for HY-80 steel. Since HY-80 
has two percent chromium it is expected to lie between pure iron and the stainless steel 
since alloy concentration 1s a major parameter at lower tempcratures affecting the ther- 
mal conductivity [Ref 12]. Thus with a high degree of confidence the thermal 


conductivity of HY-80 steel can be extrapolated to the melting point. 
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Figure Il. Thermal Conductivity versus Temperature of Iron Based Alloys 

The theoretical limit of heat capacity is well known as 3R (6 cal/mole °K) and is not 
significantly affected by alloving, except where heats of transformations are concerned. 
No good values exist for high temperature steels in general. This problem is further 
complicated by the fact that the heat of transformation is not just a function of the 
temperature but also on the material thermal history and the present rate of temperature 
change. The difference in energy stored in a martensitic steel versus a pearlitic steel is 
not negligible. Due to these factors, a precise determination of the heat capacity of 


HY-80 steel would be of linuted value. The average value of Cp for Goldak from 0 to 
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650 degree C is 536 J/kg C while a typical value from Morris is 522 J;kg C. Thus, to 
continue to use the values of Goldak for a low carbon steel is appropriate, since the data 


provided by Morris is only to 400 degrees C. 


G. CONCLUSIONS 

The use of Runge-Kutta and the explicit form of the finite difference equations is a 
powerful tool for three dimensional modeling of welding. For the unique case where the 
temperature gradient 1s controlling the accuracy this is the fastest wav to solve this 
problem. In addition, it greatly allows expanding the number of nodes in the problem 
since the computational cost is linear. Even techniques such as Gauss-Seidel for solving 
a large sparse matrix will take longer to converge the sparser the matrix is and then only 
provide an approximate solution, with no advantage gained since the gradient still limits 
tieetimic step. 

The weld model as implemented takes advantage of the lower cost of using constant 
properties were the error in doing so 1s small. In the regions of elevated temperatures, 
variable thermal properties are employed. The data available for the temperature de- 
pendence of the properties of steel is adequate for reasonable accuracy. The use of em- 
pirical weld pool results is used in lieu of simultaneously solving the weld pool dynamics. 
The later problem is compounded by the fact the altering the weld pool dynamics is the 
major reason for welding rod doping. Due to this, no advantage is to be gained by at- 
tempting to model the weld pool dynamics, since the final weld pool shape is controlled 
by the metallurgist. Thus, for the weld parameters at which it was tested, the assump- 
tions are valid and the resulting temperature profiles simulate the experimental results. 

Because of the models inherent simple structure it 1s straight foreword to alter. 
There is no mesh generator, each element can be directly altered by the user to allow 
simulating Weld flaws, torch to weld seam misalignment and other items of interest. By 
not using symmetry in developing the model, the program 1s able to study these asym- 
metrical problems. The use of explicit finite difference allowed the nodes to be described 
by three dimensional arrays, so that there 1s no nodal numbering problem which is nor- 
mally associated with using implicit techniques. This lets the user directly see what is 
happening at each node in the model, as well as the ability for simply changing its con- 
dition. This fact was exploited in the modified versions of the weld program, where the 


modifications were merely inserted and required little additional programming. 


IV. SEAM TRACKING DURING WELDING 


A. BACKGROUND 

Infrared optical svstems have been used successfully for weld seam tracking by de- 
tecting the surface thermal profiles. This is discussed in detail by Begin in (Ref. 13: pp. 
518-522] and Khan in (Ref. 14: pp. 799-805]. In both cases this was experimental work 
With no computer modeling. Since experimental results were available for comparison 


it was desired to model the occurrence of the welding arc misalignment. 


B. MODIFICATIONS TO THE WELDING MODEL 

To allow the occurrence of arc misalignment three basic steps were taken. First, a 
weld seam had to be simulated in the material. Second, the arc must reach steady state 
and then be programmed to track off to one side of the seam. And third, a run was 
conducted with a track off, with no seam present, so that the change in direction could 
be removed as a factor in the resulting thermal profile. 

The weld seam was simulated as a butt joint of two plates. The interface was a re- 
sion of zero thermal conductivity. This was created only in the fine zone region and 
was done by placing a logic statement to check to see if the region in the fine zone 
should have zero thermal conductivity or not. This check also tracked the thermal his- 
torv of the nodes adjacent to the butt joint. If the temperature of either node had ex- 
ceeded the melting temperature of the metal, the thermal conductivity was returned to 
normal, simulating that the material had fused. 

The modified programs for the model are appended with MA for misalignment. The 
FINMA subroutine 1s altered to include a logical variable called MELT which 1s used 
to determine which portions of the seam have melted. Since the seam is simulated be- 
tween nodes with x locations of 13 and 14, only the temperature history of these nodes 
are monitored. In Appendix A a detailed description of the model is given, the modifi- 
cations are to only two parts of the FINMA subroutine. This is for all internal nodes 
and for top surface nodes. In Appendix C is the listing of the modified code. Basically 
the fine node enthalpy changes are calculated normally for all nodes. The temperatures 
of the nodes adjacent to the seam are checked to see if either has melted, if they have, 
the logical variable MELT is set equal to true for that seam location. Then for those 
nodes adjacent to the parts of the seam which is not melted, the change in enthalpy is 


recalculated using a zero thermal conductivity across the seam. This approach is used 
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since it required the least modification of the code. The variable MELT is passed to the 
WELDMA main program to allow it to be saved for program restart. 

To simulate the misalignment the arc is started centered at a position of 13.5, di- 
rectly over the seam. The arc 1s started and tracks the seam until quasi-steadystate is 
reached. The arc is now shifted off to the right at .5 mm/s with a base speed of d mms. 
The surface thermal profile experiences a dramatic shift very soon after the arc is no 
longer centered on the seam. This is consistent with the results of experiments in seam 


tracking. 


C. MODEL SIMULATION RESULTS 

In Figure 12 the case of a weld misalignment is shown with a seam present. Here 
a very dramatic change Is seen, especially relative to the axis perpendicular to the weld 
axis. This change in the surface temperature closely matches the observed experimental 
result. It is important to note that the misalignment is first detected before the arc. In 
fact, due to the size of the weld pool. the weld quality may still be satisfactory when the 
misalignment is detected. In this case, there is time to effect corrective action. 

The model also shows that the surface temperature profile 1s dominated by the 
Shallow region near the surface. In these simulations, the weld seam 1s only 6.5 mulli- 
meters deep, the rest of the 20.5 millimeters is solid. This implies that even for multipass 
welds, a shallow seam is adequate for seam tracking. It also indicates that tracking the 
seam successfully says little about weld penetration. This later concern has been studied, 
and it was found that weld penetration could be correlated to weld pool diameter and 


the peak weld pool temperature. 


D. CONCLUSIONS 

The major effect noted is a rapid decrease in temperature on the side of the arc that 
the seam is on. This is due to the poor conduction across the butt joint if melting does 
not occur. This effect is important since verv primitive thermal sensors are capable of 
detecting this tvpe of change. Extensive testing with flaws, as discussed in Chapter 5, 
has shown that flaws in general always cause a rise in all temperatures in the plate due 
to the overall reduction in thermal conduction. Since the only flaws which caused a 
decrease in temperature were equivalent to localized misalignments, it can be accepted 
that a rapid decrease in temperature on one side of the arc is caused by the arc moving 
off of the seam. In addition. subsurface flaws are rarely detected prior to the arrival of 
the arc. By monitoring the temperature ahead of the arc, corrective action can then be 


correctly taken to redirect the arc toward the colder side. This is in agreement with the 
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could be used in lieu of an expensive and delicate camera system. This would provide 


a low cost and less intrusive approach. This concept is discussed further in Chapter 7. 
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V. FLAW DETECTION DURING WELDING 


A. BACKGROUND 

The prevention and detection of weld defects is the major source of cost associated 
with large scale arc welding. Due to the low metal deposition rates that are allowed 
when GMA welding the new HSLA series steels, Weld quality becomes even more crit- 
ical. To support the expert welding system, data must be acquired to interpret the me- 
aning of the various thermal signatures that any infrared imaging system would produce 
of the arc Welding process. It is desired to know the detectability of a given subsurface 
flaw for a given set of welding conditions to determine if the flaw is detectable, and if so, 
how does it change the quasi-steadystate thermal profile? 

Little experimental work has been done in flaw detection using a infrared camera 
svstem. Khan in [Ref. 14: p. 801] performed experiments in detecting surface defects and 
reported that a small A/,O, particle on the surface was detectable by a distortion of the 
temperature isotherms of the weld pool. These perturbations would start as early as 1.5 
inches before the center of the arc passed over the flaw. Since the low thermal conduc- 
tivitv impurity was on the surface, it always shows as a cold spot in the surface temper- 
ature field. 

No analvtical work has been performed in flaw detection. Current weld models rely 
heavily on symmetry, and flaws are generally asvmmetrical. Since the current model 
does not use symmetry, it 1s ideally suited for studying this type of problem. A large 
number of flaws was simulated using a modified code, which 1s designated with the suf- 
fix, LF, for Lack of Fusion. These data runs were then reduced to develop a single 


correlation which could be used to predict the detectability of subsurface flaws. 


B. MODIFICATION OF THE WELDING MODEL 

The flaw is simulated by a group of nodes which had zero thermal conductivity. 
The location of these nodes was restricted to the internal and top surface nodes of the 
fine zone for ease of programming. The basic source code 1s explained in detail in Ap- 
pendix A and the modified source code is listed in Appendix C. The flaw is rectangular 
in shape and its size and shape are specified by inputting the X and Z coordinates of the 
corners of the flaw plus its length in the Y direction. These five input variables are la- 
beled ]1, 12, K], K2, and JL. I1 and I2 are the limits of the flaw in the X direction, KI 
and K2 are the limits in the Z direction and JL is the length of the flaw in the Y 
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direction. Note, since the arc moves in the Y direction, there are no fixed Y limits. The 
relative position of the arc to the flaw is determined by a variable called LOF (Lack Of 
Fusion). LOF is calculated by equation 5.1, 


LOF = 65 — INT(VEL x TIME) (5.1) 


which is for a velocity of 4 mm’s, and a time of 10 seconds. Thus at time 10, the flaw 
1s at position 25 on the fine zone plate. For different times to quasi-steadystate and 
different weld speeds, equation 5.1 should be modified in the source code accordingly. 

The subroutine FINLF subroutine is modified as little as possible for ease of pro- 
gramming. The thermal conductivity at each node is calculated for the fine zone and 
saved in a new array called Ck. Then the thermal conductivity for those nodes presently 
at the flaw location are zeroed, using the variable LOF as the position reference. The 
location of the flaw is limited to nodes withi=3 to 25,j=3 to 25andk=1 to 6. This 
is because only the surface and center nodes were modified to allow flaws. Thus, for just 
the internal and top surface nodes, a new form of the difference equation is used, using 
a function called HK vice GK defined by equation (A.4). 


i 


CAG + LyjA) + CA(xy.z) 


De Kee elem oles) 
HK + 1.4.x v2) = cee Se Seeing oe ae 


(5.2) 
Equation (5.2) defines Hk, the major difference 1s that it uses the precalculated thermal! 
conductivities vice recalculating them as GK does. 

The main program, WELDLF, 1s modified to allow starting the program with the 
inputs as the size and location of the flaw. This program was usually run for only 5 
second simulations, with several being done in a night. A typical EXEC 1s shown in 
appendix C for performing this multiple operation using VMSCHED. WELDLEF first 
produces FILE VERIF and echoes the input and base line data for future reference. A 
second new file is created called FILE COMP, for a comparison file. Every quarter 
second of simulation an arc profile is saved. The profile is taken at the node which 1s 
three millimeters behind the arc, since experimentation showed that this would be the 
good place for monitoring the effects of flaws. This data was then processed further in 


a number of ways to determine the effect of the flaw. 


C. COMPUTER MODEL RESULTS 
In Figure 13 are contour plots showing the effect of a tvpical surface flaw. The 


torch has been traveling for 10 seconds when the flaw is suddenly introduced 5 mm 
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ahead of the arc center. The large drop in the local temperatures near the flaw is clearly 
seen from the distortion in the contours as the arc moves bv it. The deviation in tem- 
peratures 1s observable ahead of the flaw and could be used for real time process control. 
Do to an error in the program, the results that follow used a thermal conductivity only 
half of that desired for the calculations in the fine zone. Sample runs with the correct 
thermal conductivity shows that the following results remained qualitatively correct. 

The surface temperature changes from the established pattern, resulting due to a 
relatively deeply embedded non-fusion zone are seen in Figure 14. The temperature 
change from the quasi-steady levels are seen 3 mm behind the center of the heat source 
parallel to the x axis for four different flaw sizes. All flaws begin at x=11 mm andz=4 
mm and have Ax=3 mm and Az=! mm. At t=10 seconds the front surface of each 
inclusion is 9 mm ahead of the arc center. The flaw lengths in millimeters are; (a) 
Av=4, (b) Av=3, (c) Av=3 and (d) Av=1. The arc power is 2544 watts and the torch 
velocity 1s 4 millimeters per second. This simulates the response of a linear arrav of 
non-contact sensors attached to the welding torch, sensing surface temperatures behind 
the arc. 

Since these flaws are well below the surface, they are detectable by changes in the 
surface thermal profile onlv after the arc center passes over the flaw. They are then de- 
tected as a local hot spot, apparently due to the lack of conduction of heat through the 
flaw, raising the nearby interior and surface temperatures. It 1s clear from comparing 
Figure 13 and Figure 14 that unlike surface flaws, interior material defects such as slag 
inclusions or a simple lack of fusion are not easily detected prior to the arc passing over 
them. Thus while these flaws cannot be prevented by early detection, they mav be de- 
tected when the arc passes over them. The detection and marking of sites of possible 
flaws is still a useful aid in improving the weld quality. 

The temperature deviations on the surface due to near surface, embedded flaws are 
seen in Figure 15. All four flaws are located 3 mm to the left of the arc path and begin 
2 mm below the surface. Each has Ax = ! mm and Az = 2 mm. The flaw size in the 
y direction, Ay, decreases from 4 mm to 1 mm in 1 mm increments as in Figure 14. For 
all four flaw sizes, there 1s now a local hot spot between the area and the flaw and a 
larger local cold spot on the opposite side of the flaw. Temperatures between the area 
and flaw rise temporarily due to the local accumulation of thermal energy. On the op- 
posite side of the flaw the temperatures drop, since the flaw provides a constriction re- 


sistance to heat flow. The power and velocity for this example is the same as before. 
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Figure 14. Deep Flaw Modification to the Surface Temperature Profile 
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D. CORRELATION OF FLAW GEOMETRY TO SURFACE TEMPERATURE 
CHANGES 

For a given material and weld parameters it should be possible to establish a corre- 
lation between the flaw size and location to its effect on the surface temperature. For 
plain carbon steel a series of different flaws and geometries were run to develop this 
correlation. The variables that could effect the surface temperature profile are listed in 
Table 1. The material properties, arc power and velocity were held constant for each 
series of trials. The distance behind the arc was selected as a compromise on time re- 
sponse and observed change in temperatures. For the following result, a scan position 
of three millimeters behind the arc was used. The other five variables were varied one 
at atime. The resulting correlation 1s shown in equation (5.3) and is based on a total 


of 32 different combinations of flaws and geometries. 





loss ee oe) 


[om Ire 





max 
where 


Ay = AXA Y pam? 
R= S V?+Z? mm 
7,5; = The quasi-steadystate temperature at an 
x-y scan location on the surface in 
decrees C, refetenced to the 1mitiaiemiper tule. 
7 = The measured temperature at an x-v scan 
location on the surface in degrees C. 
J = 23.mm, The correlation constant for Q = 2544 watts, 


low carbon steel, v = 4 mm/s, ¥, = 3 mm. 


The left-hand term of equation (5.3) was selected since it is easy to measure. The tem- 
perature profile three millimeters behind the arc is periodically scanned. The maximum 
relative change at anv X location from the quasi-steadystate value of the temperature 
was found to be an excellent measure of existence of a subsurface flaw. The absolute 
value is for the case of shallow flaws which cause a larger local cold spot. In this case, 
the maximum temperature change that occurs 1s negative. This correlation has three 


major components: 


e The detectabilitv of a given flaw drops off as the cube of its depth in the material. 
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Table 1. SURFACE TEMPERATURE VARIABLES 


PQ | Theheatinputrate 
The material density 
The thermal heat capacity 


The depth of the flaw (measured from flaw top) 
The width of the flaw 
The length of the flaw 
The height of the flaw 


e The detectability of the flaw drops off inversely to the radial distance of the flaw 
from the arc. This is very similar to Rosenthal’s analytical solution for a point 
source. 





e The detectability is directly related to the horizontal area (XY) of the flaw. 


The surprising result was that the depth of the flaw, AZ, had little effect on the 
surface temperature profile. This is due to the dominance of the Z cubed term. As the 
flaw extends deeper into the plate, the effect of the extension drops off as the cube. 
Thus, this term can be neglected with little loss of accuracv. The above correlation 1s 
accurate Within about 20 percent for flaws ranging from about | cubic millimeter to 50 
cubic millimeters. It is theorized that this constant is accurate for similar welding con- 
ditions that have identical o and heat input rate, g/v. The difficulty in extending the 
model is that the thermal properties are verv nonlinear exactly at the point Were the best 
temperature changes are occurring. The correlation would be much more accurate if the 
temperatures were monitored far from the arc. However, the most sensitive region to 
flaws, are those areas which have the highest thermal gradients. This 1s because the 


thermal gradient acts as an amplifier of the perturbation caused by the flaw. 


FE. EXPERIMENTAL VERIFICATION 
A simple experiment was run to attempt to verifv a the predictions about the de- 


tectability of flaws. A one inch plate of HY-80 steel 7 by 12 inches was used, similar to 
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the welding model. To simulate a flaw of four square millimeters, 3/16 inch holes were 
drilled into the plate from the back side to just below the surface. The remaining depths 
tested were 0.5 and 1.5 millimeters. It was felt that these holes would be an excellent 
simulation, since the model predicted that the sensitivity would be dominated by the 
depth below the surface, and not the extent of the flaw. A GTA weld of about 2500 
Watts was passed by the flaw at about three millimeters per second. An infrared camera 
with five degree Celsius sensitivity was used to monitor the arc thermal profile as it 
passed by the holes. 

As the model predicted, the results for this case were hard to detect. Some minor 
deviation in the thermal profile was noted when the arc passed by the flaw. The shal- 
lowest flaw had the largest effect, being only a rise in temperature. The predicted drop 
in temperature was not observed, probably due to the flaw passing into the weld pool. 
The camera only produced thermal contours and like those produce by the model. these 
showed little effect from a local subsurface flaw. The temperature changes predicted by 
the model were on the order of 10 to 50 degrees. To effectively detect this, thesmmaael 
simulated a scanner at a fixed distance behind the arc. The use of an actual device per- 
forming this, instead of an infrared camera, will be necessary to reliably detect subsur- 


face flaws. 


F. CONCLUSIONS 

The model can be used to generate correlations between flaws and the surface tem- 
perature profile for any given material and welding parameters. For this particular case, 
a one cubic millimeter flaw located 4 mm down in the plate and 4 mm off of the arc 
center line would cause onlv a 0.3 percent change in the any portion of the temperature 
profile as the arc passes bv the flaw. A flaw this small and at this location would be lost 
in the noise of a real welding process. However, larger flaws would still be detected. 

The modeling of flaws during welding has provided some additional insight into the 
detectability of flaws. In general, the surface temperature profile will only detect flaws 
that are associated with the actual weld, such as slag inclusions, in and around the heat 
affected zone. This is due to its inability to detect beyond several millimeters into the 
work piece. Such detection and marking of flaws would be valuable for real time moni- 
toring of weld quality and as an aid to post weld inspection. However, due to its limited 
ability to detect flaws, it is unlikely that it could be used in lieu of current inspection 


techniques. 
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VI. THE MEASUREMENT AND CONTROL OF COOLING RATES 


Weld quality is dependant upon the cooling rate of the weld metal. The major three 
parameters of concern are; the vield strength, the nill ductility toughness and hvdrogen 
embrittlement (delaved cracking). The yield strength increases with higher cooling rates. 
The mill ductilitv toughness 1s material dependant. For example, in HY-80 higher cool- 
ing rates increase the toughness, while for HY-130, higher cooling rates reduce the 
toughness. The higher the cooling rates, the more likely hydrogen is to become trapped 
in the weld metal. increasing brittleness, which for HY-80 can appear as delaved crack- 
ing. The new high strength, low alloy steels require even tighter control of their cooling 
rates than the HY series to have satisfactory weld quality.[Refs. 11: pp. 18-21, 15] 
Mr. Morris at David Tavlor Research Laboratory requested that the computer mo- 
del be used to study cooling rates. He was particularly interested in controlling cooling 
rates during the period of arc start up and shut down. During start up and shut down 
the cooling rates are far from the steady state values and the material properties are 
usually unacceptable. Industrial practice is often to have run-on and run-off pieces of 
material adjacent to the desired weld piece which will be removed after the weld. For 
automatic welding, this is not often possible or too costly, in which case the start and 
stop areas are ground out and rewelded by hand. When welding by hand the 
arc starting 1s also verv important. The arc and puddle do not have the full pro- 
tection of the electrode coating at the instant of starting and porosity can result. 
A proven and recommended starting procedure is to strike the arc an inch or so 
ahead and then rapidly back-step to the desired starting spot. In this way, the small 
amount of metal deposited during the start will be remelted as welding progresses 
and cleansed of any gas or impurities.[Ref. 1] p.1] 

A. THE MODIFICATIONS TO THE WELD MODEL 

The original weld program did not directly measure cooling rates. In addition, the 
temperatures that cooling rates are measured are typically 30 millimeters behind the arc. 
This required that the welding model be modified. The first modification was to insert 
a subroutine which would measure and report the cooling rates. The second modifica- 
tion was to increase the size of the fine and medium zones so that most of the meas- 
urement would occur in the fine zone, which has the highest accuracy. This modification 


approximately doubled the number of nodes in the model and hence doubled the 


4] 


computational effort. The modified version is called WELDC and 1s included in Ap- 


pendix C with a more detailed explanation of the modification. 


B. CORRELATING COOLING RATES TO MATERIAL PROPERTIES 
The first question is how to measure the cooling rate. There are three commonly 


used standards, all three of which were measured: 


¢ The first is the instantaneous cooling rate at 535 degrees Celsius. This is measured 
by taking the average over the range of 550 to 520 degrees Celsius. It is considered 
the baseline cooling rate for deternuning the properties of the HY series steels. 


e The second is the cooling rate at 650 degrees Celsius. This is measured by taking 
the average over the range of 800 to 500 degrees Celsius. It is correlated to the 
effect of hydrogen cracking. 


e The third is the cooling rate at 325 degrees Celsius. This is measured by taking the 
average over the range of 400 to 250 degrees Celsius. It is correlated to the material 
strength. 


The limits for these cooling rates are based on measuring them along the arc centerline. 
Obviously, as one moves away from the center line the cooling rates will be less. It 
should be noted that even though the cooling rate at 535 degrees Celsius is centered at 
a lower temperature than the one at 650 degrees Celsius, it in general will have the higher 


cooling rate because the band it 1s the average of 1s so much smaller. 


C. THE COOLING RATES DURING STARTUP AND SHUTDOWN 

To study the cooling rates an arc power of 1950 watts and a velocity of 4 mm’s were 
selected after some experimentation. This is a fairly slow speed and low power arc, and 
hence has higher than normally desired cooling rates. Typical production weld cooling 
rates are about 25 to 20 degrees Celsius per second at 535 degrees Celsius. For the low 
power used. the cooling rates measured were typically 150 degrees Celsius per second. 
The reason for using the lower power is so the metal will have cooled to the desired 
range before leaving the medium zone of the model. Even at this lowered power it was 
necessary to double the size of the medium and fine zones to obtain this. In addition, 
the higher the power and the slower the cooling rate, the longer the model simulation is 
required to run to obtain a complete set of results. For the runs presented below, each 
was a simulation of 30 seconds of welding and cooling. Each run of 30 seconds took 
approximately six hours of CPU time to complete. The results of this numerical testing 
are intended to see how cooling rates during transients in general may be controlled. 


Due to a programmung error. the thermal conductivity in the fine zone was only one half 
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Figure 16. Baseline Velocity and Power Program 


of the required. A test run was conducted using the correct thermal conductivity and 
verified that the following results were qualitatively correct. 

The baselme power and velocity program shown in Figure 16 is a simple arc start 
and stop. The resulting cooling rates are shown in Figure 17. Ata position of zero, the 
arc 1s turned on and commences to move to the right at 4 mm/s. After a run of 80 mm, 
the arc is extinguished. At the start of the arc there are very high cooling rates. This 
is because little energy has yet diffused in the material and hence the thermal gradients 
are very steep. After the arc has traveled 15 mm, quassi-steadystate is reached and 
cooling rates are steady. When the arc shuts down, there is a rapid increase in cooling 
rates again, much for the same reason as at the start. In Figure 17 are also the cooling 
rates at 535 degrees Celsius for several distances off of the centerline. As expected, as 
you move away from the arc, the thermal gradient is less steep and the cooling rates go 
down. These lower cooling rates do not always adversely effect material properties, since 
the time above the transition temperature is less further out into the heat affected zone. 


A more graphic display of the pattern of cooling rates is given in Figure 18 where a 
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contour plot of the cooling rates around the are path is shown for 650 degrees Celsius. 


This clearly shows the unacceptable cooling rates in the vicinity of the arc stop and start. 
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Figure 18. Baseline Contour Plot of Cooling Rate: The cooling rate is measured 
Te coUrdecerc, 


D. CONTROLLING COOLING RATES BY WELD SPEED AND POWER 
CHANGES 

To control the cooling rate we can change both arc speed and arc velocity. Before 
randomly changing these parameters, the baseline profile was exanuned. It was reasoned 
that the high cooling rates at the start and stop of the welds were due to an energy de- 
ficiency. Therefore the best solution would be to augment the energy at the ends. Be- 
cause the protective shield gases are still forming at the start of the weld it was decided 
to have the arc start slowly and accelerate to the normal operating speed. To test these 
ideas, a series of four different programs were tried using different combinations of torch 
speed and arc power. 

1. Program One, Velocity Ramps 

A constant acceleration was chosen for the first try. Since the quasi-steadystate 

took about ten millimeters of arc travel to reach, this distance was selected as the ac- 
celeration range. For the baseline this distance took 2.5 seconds, using constant accel- 
eration it takes 5.0 seconds. Using this and constant arc power, the heat added over the 
first ten nullimeters would be doubled, with more of the heat being added at the start 
of the arc. This would seem to nicely make up the energy deficiency. Similar reasoning 
was used for slowing the arc down at a constant acceleration prior to extinguishing the 


arc. This programming is shown in Figure 19. 
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Figure 19. Velocity and Power Program One 


The results of program one are shown in Figure 20. As desired, the cooling 
rates at the arc stop are now near the quasi-steadystate, the 650 degree Celsius cooling 
rate having been corrected the most. At first glance, it would appear that the drop in 
cooling rates at the very start and the trailing edge of the arc is a problem. This in in 
fact the correct response. Looking at the cooling rates off of the arc centerline at 535 
degrees Celcius can help clarify this. The centerline cooling rate has almost dropped to 
the cooling rate five millimeters off of the centerline at the start. This occurs at a dis- 
tance of minus five nullimeters behind the torch starting position. Because the torch has 
never crossed this point, it is supposed to have a cooling rate similar to the heat affected 
zone. In Figure 21 this is clearly seen. The cooling rate contours with the programming 
are almost symmetrical around the torch path. At the start, the cooling rates at a fixed 
radial distance from the arc start location are the same. 

Having the arc slow down at the end of the run to add the extra energv did not 
have the desired effect. The cooling rate at the point of arc slowdown dramatically 


dipped and the cooling rate at the point of arc off still steeply climbed and a slightly 
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Figure 21. Program One Contour Plot of Cooling Rate: The cooling rate is mea- 
sured at Op deg ne: 


higher cooling rate were obtained. The onlv benefit gained is that the cooling rates at 
stop also now show a radial symmetry. Reconsidering the problem, it is noted that the 
energy deficiency 1s not in the region near a quasi-steadystate profile, but effectively at 
the point the arc stops. Thus to correct this, maybe power, and not velocity program- 
mung should be used. 
2. Program Two, Velocity Ramp on Start, Power Ramp on Stop 

The velocity and power program two Is shown in Figure 22. Here, once the arc 
has traveled the desired distance, the torch stops and the arc is left on. The arc power 
is linearly reduced to zero over a 2.5 second interval to make up the energy deficiency. 
Only one half of the estimated energy deficiency is made up in an attempt to control the 
cooling rate dip prior to turning off the arc. In (b) of Figure 22 it can be seen that the 
results are nuxed. The dip in cooling rate prior to securing the arc is mininuzed and the 
650 degree cooling rate is noticeably improved. However, the 535 degree Celsius cooling 
rate has slightly increased. This later effect is due to the way these cooling rates are 
measured; the 650 degree Celsius cooling rate is averaged over a 300 degree Celsius band, 
While the 535 degree Celsius cooling rate is averaged over a 30 degree band. This indi- 
cates that the ability to control a cooling rate depends on what temperatures it is taken 
Over 
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3. Program Three, Velocity Ramp and Holding Power 
To further attempt control over the 650 degree Celsius cooling rate, the program 
shown in Figure 23 was tried. This starts as the others, but the power is held constant 
after the torch stops for an additional 2.5 seconds, vice ramping down as before. This 
would completely make up the estimated energy deficiency. As can be seen in part (b) 
of Figure 23 the 650 Celsius degree cooling rate was lowered even closer to quasi- 
steadyvstate when the arc is secured. As expected, the dip in cooling rate prior to stop- 
ping the arc has increased. The maximum 535 degree Celsius cooling rate is identical to 
the baseline. 
4. Program Four, Constant Power to Velocity Ratio 
Since in the quasi-steadystate, the cooling rate is directly rated to the ratio of 
arc power to torch velocity, it was conjectured that keeping this ratio constant while 
securing the arc would result in a constant cooling rate. This programming is shown in 
Figure 24. As can be seen in (b) of Figure 24 this is not the case. This was in fact the 
worst of the four programs with the largest excursions in both the 535 and 659 degree 


Celsius cooling rates. 


E. CONCLUSIONS 

It is possible to change the cooling rate during transients such as arc startup and 
shutdown. The best results Were obtained for controlling the cooling rate when starting 
the arc. Since there 1s no thermal history affecting the diffusion of heat, significant 
changes can be affected as desired. The control of cooling rate on arc stop is more dif- 
ficult. In this case there is a thermal history which will effect the cooling rates. The 
major noted problem, is that trving to make up an estimated energy deficiency caused a 
dip in the cooling rate prior to stopping the torch. This is because adding more energy 
at the end of the weld increases the weld pool size. Thus the region of the arc next to 
this enlarged pool has less low temperature material than the quasi-steadvstate portion 
of the weld bead in which to diffuse heat. This results in lower cooling rates. Of the 
tests performed, program two seemed to be the best compromise. 

Another important result is the idea of energy deficiency. By use of a simple energy 
balance and knowing the distance to quasi-steadystate it was possible to select a simple 
velocity and a power programming to keep the cooling rates near their quasi-steadyvstate 
values. This is important for the distance to quasi-steady state can be easily determined. 
Using this to develop a velocitv and power program provides two major advantages. 


The first 1s that the run-on and run-off pieces can be reduced or eliminated. The second 
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1s that expert welders designed to control the cooling rate can use preprogrammed 
Startup and shutdown algorithms, since during these transients anv system will have in- 


adequate feedback for effective control. 
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VII. THE INFRARED DIODE DETECTOR 


A. INTRODUCTION 

For many vears Welding was a skilled trade that was at times more art than science. 
As the study of the metallurgy of welding progressed, new welding techniques were in- 
vented. With the advent of modern robotic controls, welding tasks began to become 
automated. But, due to an incomplete understanding of all that occurs in the weld pool 
during the welding process, designing a closed loop control system of an automated 
Welder is not a simple task. 

Current automate Welding svstems attempt to eliminate as many variables as possi- 


ble in the welding process, to simplify the task, and ensure better production rates. But 


.. these svstems are incapable of correcting for perturbations that arise during the 

welding process. Control of the welding process requires the identification and 
monitoring of perturbations in a wide variety of parameters. The varied nature of 
these parameters and the large number of variables involved have thwarted previous 
attempts at closed loop control of the welding process.[Ref. 16: p. 799] 


But even when care has been taken to remove the obvious variables, ensuring that there 
is a near perfect fit with no gaps between the pieces to be welded, for example. problems 
still arise. 

The Edison Welding Institute reported in Mav, 1986, on a research project dealing 


with robotic vision in automated welding systems that 


..accumulated variations in the motion and plate movement due to warpage and 
residual stress may cause the robot to miss the intended weld joint. The solution 1s 
to give the robot an eve and other sensors which “see” and “feel” when a correct weld 
is being made and can adjust itself when deviations occur.[Ref. 17] 


This has been performed successfully by several researchers. Khan states that 


..attempts to adapt intelligent vision systems to seam tracking and weld puddle 
contro] have been made. Most attempts of welding control have been adapted from 
computer-vision-based systems...recently infrared thermography has been used to 
monitor cooling rates in Welds as a possible means of on-line-control of heat input. 
Although the above research has lead to increased knowledge of arc welding physics, 
an acceptable system for in-processing welding and quality control has not been 
developed. The development of improved sensors is required....[Ref. 16: p. 799] 


A review of the literature found three different infrared thermography systems in 


use. These are the infrared scanning camera, the fiber optic spot sensor, and the solid 
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state video camera. Of these three, R. A. Morris, from David Taylor Research, reports 
that the solid state video camera is preferable to the nitrogen cooled scanning camera 
due to “...ease of data acquisition and reduced camera alignment problems.’ [Ref. 15] In 
addition the camera is less expensive to purchase and to operate. Given the require- 
ments for a reliable infrared sensor, using infrared diodes was suggested as a low cost 
alternative. Because thev are equivalent to a non-contact thermol couple, the diodes 
would provide much higher data acquisition rates, one area in which all present sensors 
are only marginal. This chapter explores the feasibility of using infrared diodes to 


monitor the weld process. 


B. CHARACTERISTICS OF INDIUM ARSENIDE DIODES 

The temperature band of interest for monitoring Weld cooling rates 1s centered about 
800 degrees Kelvin. Temperatures above this, at about 1300 degrees Kelvin, are useful 
for monitoring the heat affected zone to detect flaws and welding arc misalignment. The 
peak weld pool temperatures are around 2400 degrees Kelvin. A search of technical lit- 
erature on infrared diodes located the Indium Arsenide series as having the best spectral 
response at room temperature. The response 1s shown in Figure 25, the peak radiative 
temperature corresponding to 2.0 microns and 4.0 microns are 2762 degrees Kelvin and 
1381 degrees Kelvin, respectively. However at these low temperatures, the curves are in 
general verv flat. That is, at low temperatures, the energy 1s spread out over a large 
frequency band. Thus, based on spectral response, the indium arsenide diodes are an 
acceptable sensor for detecting the temperature range of interest. 

The indium arsenide diodes have a very small junction resistance, typically about 50 
ohms. Thus, thev are onlv operated in the short circuit mode. In the short circuit mode 
thev provide linear current output versus light input for up to 9 order of magnitude. In 
addition, in the short circuit mode, the response of the diode 1s essentially independent 
of the diode temperature, removing the need for temperature compensation of the de- 
tector. Ref. 18] 

The indium arsenide diode purchased for testing was an ORIEL 71110, with a pur- 
chase price of $275.00. This diode has an active area of .79 mm? For testing it was 
mounted in ORIEL 7192 detector module, which is for unbiased operation. The detec- 
tor provides a reverse current protection diode, since that maximum safe reverse voltage 


On an indium arsenide diode is one volt. 
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Figure 25. The Spectral Response of Indium Arsenide Detectors 


CO THE DIODE DETECTOR ASS EME 

Other than the window of the diode, there is little control over the directionality. 
While advances have been made in infrared optics, a review of the literature showed that 
use of simple light pipes would be more than adequate for the type of focusing desired. 
Based on the results of computer modeling of flaws and misalignments, it was deter- 
mined that a spatial averaging of the temperatures over an area of several square milh- 
nieters was a desirable feature of the sensor. Light pipes were sclected sinee their use 
“instead of a mirror and lens opties results in increased economy and simplicity of the 
apparatus. ”"[Ref. 19] 

The design of the hight pipe was based on three considerations. To keep the diode 
out of the wav of the welding torch. To minimize the signal loss due to the length of the 
light pipe. And to have an acceptable spatial resolution. After several trial light pipes, 
the one shown in Figure 27 was selected. A series of initial calibration runs were per- 
formed to find the best cireuit. A zero bias, zero offset as shown in Figure 26 was 1ni- 
tially used, with a gain resistance of Rg = 2K. Prior to calibration the offset output 
voltage was adjusted to a minimum by the trim pot, Ro. Since this particular cireuit 
used an inexpensive 704 operational amplifier, the offset would drift as the circuit 
warmed up. Use of higher quality components would eliminate this problem. The drift 


in offset docs not effect the measured output of the diode since it is always possible to 
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Figure 26. A Typical Zero Bias, Zero Offset Sensing Circuit: The circuit 1s 


shown using a 741 operational amplifier. 
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Figure 27. The Infrared Detector with Light Pipe: The stainless steel light pipe 
was press fitted into an aluminum holder, with the diode resting in its 
own mounting disk. Both are held into the Oriel 7192 detector with set 
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Figure 28. Sensitivity of a Light Pipe Detector to Distance: The calibration was 


done using the zero bias, zero offset sensing circuit. 


measure the offset and subtract it awav. What it does do, is limit the minimum tem- 
perature which can be detected, since the diode signal must be greater than the offset to 
be detected: 

In Figure 28 is shown how the light pipe minimizes the sensitivity of the detector 
to minor changes in distance from the heat source. This 1s important since it makes 
alignment less critical and removes one more variable from the detector calibration. As 
expected, the light pipe did provide adequate focus so that little signal strength was lost 
for minor changes in distance from the source. A calibration curve is shown in 
Figure 29. The non-linearity does make the detector less suitable for measuring of 
cooling rates, however it is still suitable for seam tracking and flaw detection. 

When the diode was zero biased, it was verified that the detector was insensitive to 
temperature changes. This was important, for in the biased condition the detector will 
rapidly lose signal strength as it heats up, due to its change in internal resistance. Test- 
ing showed that the diode would heat at as much as 0.5 degrees Celsius per minute. This 
heating was due to the radiant energy being absorbed by the diode. When the diode was 


tested in the open circuit mode, versus the short circuit mode, raising the temperature 
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Figure 29. Typical Calibration Curve for an Indium Arsenide Diode: The cali- 
bration was done using the zero bias, zero offset sensing circuit. The 


diode dark voltage was five millivolts. 


from 10 to 30 degrees C resulted in a 75 percent loss of signal. Thus, all further exper- 


Imentation was done in the zero biased mode, that is short circuited. 


D. CONCEPTUAL IMPLEMENTATION OF INDIUM ARSENIDE DIODES 

The use of a simplified infrared vision scheme using lead sulfide photo conductance 
detectors was developed by Begin [Ref. 13: pp. 520-522]. This system was designed 
specifically for seam tracking and had a three dimensional capability. This allowed him 
to not only control the torch position relative to the seam, but also the torch height to 
control the arc length. The detection system he developed used a fiber optic cable to 
transmit the infrared energy from the work piece to his array of lead sulfide detectors. 
Though not mentioned in his article, this appears to be due to the bulkiness of the lead 
sulfide array and the need to keep these conductance detectors at a constant temperature 
to maintain calibration. 

The indium arsenide diodes can directly replace the more complicated and expensive 
fiber optic cable and lead sulfide detectors. Since the diodes are only 1/8 inch in diam- 


eter, they can be directly mounted in to a light pipe assembly, such as shown in 
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Figure 30. Proposed Configuration of a Light Pipe Array 
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Figure 31. Proposed Mounting of a Light Pipe Detector Array 
Figure 30, and still be in a package small enough to mount onto the weld torch. A 


possible installation of such an array is shown in Figure 3i. Further more, it would be 
possible to mount individual diodes at exact positions of interest. This is useful, since 
no one location is ideally suited for taking every type of measurement. A list of possible 


locations are: 


60 


¢ The side of the arc: Looking for subsurface defects. 
e The leading transverse profile: Looking for torch to seam misalignment. 
e An angled leading transverse profile: Looking for torch height from work piece. 


e Behind the arc: Measuring the cooling rate for material property control. 


The spatial resolution of the detectors is a matter of some concern. Begin felt that 
the resolution was dependant only on the detection area of a given detector[Ref. 13: p. 
S21). If a detector had a spatial resolution of 0.2 mm, then that was his spatial sensi- 
tivity. What he neglected to note is that the steep temperature gradients act as an am- 
plifier. For example, the temperature gradient in the vicinity of the arc is of the order 
of 1000 degrees per inch. If the detector has a temperature resolution of 10 degrees, than 
the linear resolution would be 0.01 inch per inch of detector. So that a detector of 0.2 
inches diameter would be able to detect an arc movement of 0.002 inches. Taking into 
account this amplification effect of the temperature gradient, it 1s possible to use fewer 
detectors with wider apertures. In addition, overlapping of detector windows, which 1s 
desirable, is now feasible. The use of the non-focusing light pipes is ideal for this tvpe 


of coverage. 


E. OPERATIONAL TEST OF THE INDIUM ARSENIDE DIODE 

After calibration and sensitivity checks were completed the diode was taken to the 
welding laboratory to check its response in the presence of torch arc. An initial sensi- 
tivity check showed that the diode was not adversely affected by the torch arc. It also 
showed that the gain was too low for the temperature range of interest. 

The diode circuit gain was increased by a factor of five, and the diode was recali- 
brated. It was mounted on a traverse to allow scanning a stationary arc thermal profile. 
When the circuit was energized. the offset was noted to have drifted since calibration to 
0.06 volts from 0.01 volts. The torch arc was started with the diode directed two inches 
away. Upon attempting to read the diode, the diode had failed in the shorted condition. 


The cause of the diode failure has not been determined. 


F. CONCLUSIONS 

The use of the indium arsenide diodes or lead sulfide conductors is a suitable ap- 
proach for a production environment infrared detector. While cameras and higher re- 
solution vision schemes are useful for research in the laboratory, onlv that information 
required to control the welding process should be gathered for actual welding fabri- 


cation. Due to their small size, light weight, minimal temperature sensitivitv and high 
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speed performance the diodes are more than adequate to gather the required thermal 
information. However, if it is determined that higher resolution is required, then the use 
of large numbers of these detectors would result in a rapid loss in any savings, and svs- 
tems such as solid state video cameras should be considered. Also, further testing will 


be required to determine the suitability of the diodes for the production environment. 
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VII. RECOMMENDATIONS 


There are many areas in which this study can be continued. The model has proved 
itself adaptable to studying a wide variety of Welding problems and can provide results 
with a reasonable computational cost. The indium arsenide diode has several features 
which make it a suitable detector in the production environment. Based on these results 
the following recommendations are made: 

¢ Perform additional experimental verification of the model results. 
e Verify the spatial dependence of the welding model. 


e Investigate the effects of power level, torch velocity and material properties on flaw 
detection. 


e Rewrite the cooling rate version of the model to take advantage of symmetry. 
e Determine the transient distances for arc startup and shutdown. 


¢ Develop a simplified vision system using indium arsenide diodes. 


Some of the results of modeling have been verified by previous experimental work. 
This includes misalignments and surface flaws. The testing that is recommended. is the 
studying of subsurface flaws. The model predicted that the detectability of the flaw is 
dominated by its depth below the surface and it honzontal area. Since detectability is 
not significantly affected by the extension of the flaw into the plate, it should be possible 
to drill holes from the opposite side of the plate to various depths below the welding 
surface, and then these holes would approximate the flaws simulated in the computer 
model. The two major results desired would be to confirm that shallow flaws cause both 
a minor temperature rise and a larger temperature drop, and that deep flaws cause only 
a local temperature rise. To measure these changes, a scanning type detector vice a in- 
frared camera will be required. 

Experimental verification of the predictions on cooling rates could be done two 
ways. The first by use of an infrared camera to monitor temperatures and cooling rates 
during a weld sequence. The second is to section the material after the weld and analyze 
the heat affected zone’s metallurgy. The last would be the more time consuming meas- 
urement, but would provide a detailed picture of the thermal history of the weld and heat 
affected zone. Using a simple power and velocity program scheme, as recommended by 
a simulation, a test weld could be performed and analyzed to determine the actual ef- 


fectiveness of controlling the cooling rate. 
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The spatial dependence of the model is still of some concern. The selected one mil- 
limeter spacing was based on the experience of other researchers and simple two di- 
mensional trials. When analyzing flaws, some boundary instabilities were observed, so 
that minor damped temperature oscillations occurred at the 1mmediate location of the 
flaw. To establish this as the cause of the oscillation two choices are available. The 
boundary can be “softened” by increasing the flaws thermal conductivity from zero to 
an intermediate value, say 10. w/m°K The second is to change the grid size. Normally 
grid sizes are halved when checking spatial dependence. To do this would require eight 
times as many nodes, and due to the time step limitation, taking one fourth the time 
step. This would cause the model to run 32 times longer! Based on this, it is recom- 
mended that the grid size be reduced by a larger fraction, say 3/4. In this case the model 
would onlv run 4.2 times longer for the same simulation. 

The detectability of flaws was limited to the study of the effect of flaw location and 
geometry. The effect of power was studied over a small range, and in general no corre- 
lation was found between power changes and the relative changes of temperature from 
quasi-steadvstate. A further study of the effects of power, as well as torch velocity and 
material properties should be pursued using the model. Based on this research it should 
be possible to generate a more general correlation between a given flaw and its detecta- 
bility including these variables. 

The program used to simulate weld cooling rate, WELDC, does not invoke sym- 
metry around the arc path axis to halve the number of nodes. Because the problem be- 
ing study was symmetrical, this could have been done. While this is a significant rewrite 
of the code, it may be required if more realistic welding conditions of higher arc power 
and welding speeds are desired. These later will require extending at least the medium 
grid nodes behind the arc even further, lengthening a program which presently requires 
about one hour of CPU for everv five seconds of simulation. By using symmetry, the 
computational cost could be halved. 

One of the interesting results of the cooling rates experiment was the importance of 
the transient length during startup and shutdown. By simply knowing this distance for 
a given set of welding parameters it was possible to develop a program to control the 
cooling rates. Being able to predict these lengths would then allow programming for any 
set of welding conditions. It may be possible to determine these distances analytically, 
and if not, repeated computer experimentation could provide a basis for prediction. 

Due to the speed of the model using Runge-Kutta, 1t may now be possible to expand 


the code to include weld pool simulation. Normally, this is not possible because the 
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dynamics of the weld pool require more nodes than are available in a model. The sim- 
plest case would be the development of a stationary Weld pool. This would allow elimi- 
nating the coarse zone and shrinking the medium Zone to reduce the number of nodes. 
It would then be possible to increase the number of nodes in the weld pool and still have 
a system that could be solved in a reasonable amount of time. This 1s particularly true 
When using a mini- VAX, which could be used, since the explicit technique uses verv little 
memory compared to an implicit technique. Thus a problem of interest could be set into 
a dedicated mini-VAX and left to run for a day, at essentially no cost to the user. 

The final recommendation is that a simplified vision system be prototyped using in- 
dium arsenide diodes. To successfully use infrared techniques in production welding re- 
quires that the detectors be as simple to use and rugged as possible. These diodes are 
relatively inexpensive and can be easily mounted with the welding torch. Thev require 
no special support equipment, and with a simple zero bias detection circuit can be di- 
rectlv fed into any analog-to-digital converter for use by an expert welding svstem. 


Further testing is required to determine their suitability for the production environment. 
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APPENDIX A. DETAILED EXPLANATION OF THE WELDING 
COMPUTER MODEL 


A. EXPLANATION OF THE SOURCE CODE 

The program has four major parts. The first is the main program called WELD. 
This controls input and output, adds the energy from the arc, uses Runge-Kutta to solve 
the semi-discrete difference equations and repositions the grids relative to one another. 
The other three are the subroutines that contain the finite difference equations for each 
grid zone; the fine, medium and coarse. These subroutines are used by the Runge-Kutta 
solver to find the temperatures at the next time step. A listing of the the variables used 
and their definitions 1s found in Table 2 below. The only inconsistency is in the use of 
the arravs in the subroutines. In the subroutines the input arrays are referred to as A, 
B or C instead of AIN, BIN or CIN. The discussion that follows refers to the source 


listings. 


Table 2. VARIABLES USED AND THEIR DEFINITIONS 


re a Definition 


[AGI 36) Array of coarse grid temperatures in degrees Kelvin 


AIN(21.36) Array of coarse temperatures for input of the next Runge-Kutta 
step 

AOUT(21,36) | Array of coarse temperature changes on output of the Runge-Kutta 
step 


ASU | ASUM(21,36) | 36) | Arrav of coarse temperature change sums used bv Runge-Kutta 


eG Re) eee Array of medium grid temperatures in degrees Kelvin 


BES 2727) Array of medium temperatures for input of the next Runge-Kutta 
step 

BOWM27 27) Array of medium temperature changes on output of the Runge- 
Kutta step 


BSU M2727} Array of medium temperature change sums used by Runge-Kutta 
CO) Array of fine grid enthalpy in Joules per cubic meter 
GIN@7Z Array of fine temperatures for input of the next Runge-Kutta step 


COUT(27.27) | Array of fine enthalpy changes on output of the Runge-Kutta step 
CSUM(27,27) | Array of fine enthalpy change sums used bv Runge-Kutta 


BI(4) The Biot numbers for different grids 
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Lsed to count the number of times the fine grid has shifted, when 
it reaches three, the medium grid shifts and it is set to zero. 


The fixed thermal conductivity used by the medium grid in W mk 


Counter of which time step is being solved 


Counter of total number of steps taken, will be different from M 
—— if problem has been restarted. 


XB | Pointer which locates the medium grid relative to the coarse grid 


Pointer which locates the fine grid relative to the medium grid 
NDIV The number of divisions (or time steps) to be taken 


Ca Compared with time to decide when to output data, is then in- 
creased bv 0.5 seconds 
The amount of enthalpy change of a given control volume for each 
time step in Joules per cubic meter 


QDOT The rate of energy input from the arc in Watts 


Lsed to compare the distance the arc has traveled, when they are 
equal the fine grid shifts and it is increased by three units 




















The spatial weighting factor for arc density, varies with arc position 


The voltage of the arc in volts 


AARC Position of the arc on the fine grid along X (1) direction in milli- 
meters 


MEL The velocity of the arc in the X (1) direction in mm s 


YARC Position of the arc on the fine grid along Y (j) direction in milli- 
meters 


Bee L The velocity of the arc in the Y (j) direction in mm's 











1, Main Program Weld Fortran 
The main program, WELD, is divided into five sections, plus three functions. 


These are input and setup, adding heat from the arc, the Runge-Kutta solver, running 
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output, and final output plus setup for a restart. The three functions are for determining 
the thermal properties of the fine grid. These are thermal conductivity as a function of 
temperature, and converting to and from temperature and enthalpy. 
a. Input and Setup of a Problem 

The first section of the main program deals with input for setting up the 
problem. Once the initial validation was complete, the program was run as a batch job. 
Thus prior to each run the program 1s recompiled with the appropriate fixed conditions. 
For each run the length of the simulation 1s specified using FINI. If a new data run was 
commencing, the GOTO 100 hne is commented out and the initial conditions are set and 
the output files are opened. If this is a continuation of a previous run, then the program 
jumps to statement 100 where the data from the previous run is loaded and the output 
files are prepared. In general, the continuous output files SURF and CUT are reposi- 
tioned to their ends to allow appending additional data. The initialization of parameters 
is skipped by reentering at statement 200. Then the material properties are defined and 
the program enters the main processing loop. 

b. Adding Heat from the Arc 

The first step of the main processing loop calculates the energy added by the 
arc. It first calculates a volume weighting factor, SUM based on the arc’s present X- 
Y location. This is necessary because the amount of energy added to a nodal point 1s 
based on its distance from the arc center, which changes each time step. Thus, to keep 
the energy added each time step the same, the spatial distribution 1s normalized by the 
use of SUM. SUM and the actual energy input for this time step, Q, to determine the 
change in enthalpy of the nodes under the arc in the fine grid. Obtaining the position 
of the arc relative to the fine grid is done using the variables NARC and YARC, where 


NARC is usually constant at 14 and YARC 1s as defined 1n equation A.1. 
YARC = VEL x TIME + 73 —-9NB—3NC (An 


YARC uses the two pointers, NB and NC, to calculate the relative position of the arc 
on the fine grid from the fixed reference of the coarse grid. The value of 73 is based on 
an initial NB = 3, NC = 10 and the arc at position 16 on the fine grid. YARC is in 
millimeters. 

The Gaussian power distribution of the arc is done using equation A.2, 
which is shown for the case of a hemisphere of radius four millimeters, that is where 


a=b=c=4 mm. 
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i Oe Fae as 
O(xy.2) = Ae: Be a b2 Bs co 
= C, exp[ — 0.10625(x? +? + 2’)] 


The value of the change in enthalpy at each node 1s approximated by taking the value 
at the nodal point and assuming it 1s constant over the control volume of the node. To 
insure a proper energy balance, a weighting factor is calculated by summing the values 
Simeerat each node. SUM is then used to normalize C,, where C, = Q‘SUM. Hence at 
each time step, the same total change in enthalpy occurs. When double elipsoids are 
used, that is the power distribution is not svmetrical about the center of the arc, then it 
1s necessary to use the first form of the Gaussian power distribution. 
c. Runge-Autta Solver 

The Runge-Kutta solver requires temperatures, so the first step 1s to convert 
the fine grid enthalpy to temperature. The function T(h) performs this conversion and 
places the temperatures into the array CIN. The basic form of fourth order Runge- 


Kutta is shown in equation A.3, 


KI" = f(T") 
rzn 
ke = fi Ts ai Al ) 
rey 
ee ae SS (4.3) 
KI" = f(T" + K3") 
piel pty (K1" + 2(K2" + K3") + K4") 
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where f(7) is the semi-discrete form of the difference equation. By combining terms be- 
tween steps, it 1s possible to use only three arrays per grid, vice the five listed above. 
The input to f(7) is the IN arrav and the K is the OUT array. After each step, a new 
IN array is calculated by adding the OUT array values, and the OUT array values are 
added to the SUM array. For the fine grid, the COUT array is in enthalpy, so that the 
new CIN values are found using the function T(H). After the fourth Runge-Kutta step, 
the A, B and C arravs are updated to the next time step by adding one sixth of the 
poe VI, BSUM and CSUM arrays, respectively. 
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d. Aloving the grids 

The variable STEP is used to decide when to move the fine grid. If the 
distance traveled equals or exceeds STEP then it is time to move the fine grid. In this 
case STEP 1s increased by three mullimeters and the counter for moving the medium grid, 
BSTEP, is incremented one. Moving the fine grid 1s done in three steps. First the 
trailing edge of the grid will be now in the medium grid. These are all fine nodes with 
ajindex of 1, 2 or 3. There is one medium grid node for each 27 fine grid nodes (18 fine 
nodes on the surface). The enthalpy for the fine grid nodes are averaged together and 
converted to a temperature which 1s assigned to the equivalent medium grid node. Care 
is taken for the surface nodes, those with a k index of 1, because they only have half the 
volume of a regular node. The second step is to shift all the values of the fine grid. 
Those nodes with a j index of 1, 2 or 3 have already been accounted for, so the whole 
fine arrav 1s shifted -3 in the] index. The last step 1s to assign values to the leading edge 
nodes, those with a } index of 25, 26 or 27. They are all assigned an enthalpy equivalent 
to the temperature of the medium grid node they will now be occupying. This conver- 
sion 1s done using the function H(T). Upon completion, the pointer NC 1s incremented 
One to indicate the new location of the fine grid. 

In the case that BSTEP is now equal to three, it 1s time to shift the medium 
grid in the coarse grid. BSTEP is set equal to zero and the erid 1s again shifted in three 
steps. The first step is to average the trailing nodal points and assign them to the 
equivalent coarse grid nodal point. There are 90 medium nodal points for each coarse 
nodal point. Those nodes with a j index of 1, 2 or 3 are averaged together, taking care 
to note that the surface nodes, those with a k index of 1 or 10 have only half the volume 
of a regular node. The second step is to shift all of the nodes -3 in the } index direction, 
since those nodes with a j index of 1, 2 or 3 are already accounted for. The last step 1s 
to assign the temperature from the equivalent coarse node to the leading edge nodes, 
those with a j index of 25, 26 or 27. Upon completion, the pointer NB 1s incremented 
One to indicate the new location of the medium grid and NC 1s reset to 10 to indicate the 
new location of the fine grid. 

e. Running Output 

During the running of the program, snapshots are taken of the temperature 
every half second and saved. The data is saved in two files, SURF and CUT, which are 
then used to generate three dimensional plots and contour plots of the temperatures. 
The data saved are the surface temperatures of the fine grid, the surface temperature of 
the medium grid and the X-Z profile of the fine grid directly under the arc. For the fine 
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grid output the function 1(H) is used to convert the enthalpy to temperature. For the 
medium grid. the nodal points occupied by the fine grid do not normally have a value 
assigned. For continuity, the enthalpy of the fine grid nodes are averaged together, 
converted to temperature and assigned to their equivalent medium grid node. This is 
identical to the process used by the grid shifting section of the program. 

f. Final Output and Setup for Restart 

Lpon completion of the program the time and the final temperatures of the 

fine, medium and coarse grid are placed in a file called FINAL. The regions in the me- 
dium and coarse grid which are occupied by another grid are filled in by averaging the 
values from the finer occupying grid, as is done during grid shifting. This file is used to 
generate three dimensional and contour plots of the temperature profiles. A second file, 
RESTAR, is used for restarting the problem. It saves the parameters necessary to re- 
start the problem as well as the three arravs; A, B and C. 

g. Functions 

ite wei erognam wses three function, bh, | and fH. Whe first is Fh, 
Which models thermal conductivity as a piece-wise linear function of temperature, which 
is discontinuous onlv at melting. The second is T, which models the temperature as a 
piece-wise continuous linear function of enthalpy. The third is H, which models the 
enthalpy as a piece-wise continuous linear function of temperature. The values are the 
same as used bv Goldak [Ref: 5: pp. 587-600]. 
2. Subroutine FIN Fortran 

This subroutine contains the semi-discrete form of the explicit finite difference 
equation for the fine grid. The parameters for this grid are listed in Table 3. The sub- 
routine 1s broken into nine sections, each dealing with a different tvpe of node. These 
nodes are shown in Figure 32, where the letters in the figure correspond to the following 
sections. The program makes use of a function GK, which is shown in equation A.4, 
to improve the readabilitv of the code. The function GK takes the harmonic mean of 
the thermal conductivities and multiplies it with the temperature difference for each dif- 
ference pair. When this is done for all adjacent nodes, the change in enthalpy at that 
node is known. The enthalpy form of the equation takes into account the variation of 


the heat capacity with temperature. 


BOS Oey, 


Ghai l,)— FK(T,) + FK(T,) 


(4.4) 
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Table 3. PARAMETERS DEFINING THE FINE GRID 


Dimensions: 27mm x 27mm x 7.5mm 
Node Volume: Imm x Imm x Imm 
Atlhave (2723) 


The following constants are used: 
FO(3) = 1.0x10° At s/m? 

B1(3) =9.399x107° At Jim? KO 
BI(4) = 5.0x10° At J/meK 


FKB = 26.5 W/mK 
Kb = 53.0 WimK 
Cp = 4.0x10° Jim? K 
h = 25.0 Wim?K 
=) oD 
og =5.67x10° W/m? K" 





The function FK(T) returns the thermal conductivity for the input temperature. 

Three of the constants in Table 3 require a little explanation. The first, FO(3), 
is like the Fourier number, and is defined by FO(3) = Ar/Ax?.. The second 1s an equiv- 
alent Biot number for radiation from the surfaces. This is defined by BI(3) = ecAr/Az. 
The last is the equivalent Biot number for convection from the surfaces. This is defined 
by BI(4) = AAtj/Az. Note that the Biot numbers are with respect to Az, not Ax and that 
Az is one half of Ax for a surface node. These two equivalent Biot numbers are the 
product of the normal Fourier and Biot numbers. 

Interfacing between grids requires care. As discussed in Chapter 2 the nodal 
spacing controls the difference equation. The distance to a medium node 1s twice that 
of normal fine node spacing. Because the change in enthalpy at a node 1s related to 
1/Ax this reduces the effect by one half. This is why FKB is one half of Kb. The second 
difficulty in interfacing is keeping track of the association of fine nodal points to adja- 
cent medium grid nodal points. In the subroutine all medium array indexes are com- 
bined with a B to indicate they are B array subscripts. 

a. The Interior Nodes 
The equation for the interior node is the simplest. This equation 1s valid 


over the following range: 1 = 2 to 26,) = 2to 26 and k = 2 to 7. 
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Figure 32. The Fine Grid Node Types: The letters coresspond to the sections of 


the subroutine FIN. 


b. The Top Surface Nodes 
The difference equation for the top surface nodes is similar to an interior 
node but it has a convection term and a radiation term. Because Az is one half of Ax, 
the difference of the k index has twice the effect of the other diflerences. This equation 
is valid over the following range: 1 = 2 to 26,} = 2to 26and k=1. 
c. The Bottom Surface Nodes (Medinm Interface) 
The difference between this equation and the interior node equation js that 
the +k index term is replaced by a difference term with the medium grid. [KB takes 


care of the difference in node spacing between the fine and medium zone. This equation 
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is valid over the following range: 1 = 2 to 26,] = 2 to 26,k = 8, IB = 10 to 18, JB = 
NC to NC +8, KB = 4. 
d. The End Faces (Medium Interface) 

The end faces are similar to the bottom nodes. Depending on the face, ei- 
ther the +] or -j index term has been replaced by a medium grid difference term. As 
before, care was taken in matching the fine to the medium nodal points. These 
equations are valid over the range: 1 = 2 to 26,)] = 1 or 27, k = 2to 7, IB = 10 to 
18, JB = NC-1 or NC+9, KB = 1 to 3. 

e. The Side Faces (Medium Interface) 

The side faces are similar to the end faces, with either the +1 or -i index 
term being replace by a medium grid difference term. These equations are valid over the 
range: 1 = 1 or 27,] = 2 to 26,k = 2 to 7, IB = 9 or 19, JB = NC to N@ ema 
=F iwc Se 

f. The Top and Bottom End Edges (Medium Interface) 

This is a combination of dealing with grid interfaces and surfaces. The top 
equation has either the +] or -j] index term replaced by a medium grid difference term 
and the -k index term 1s replaced by the two surface terms for convection and radiation. 
The bottom equation has either the +] or -] index term and the +k index term both 
replaced by a medium grid difference term. These equations are valid over the range: 1 
= 2 to 26,j = 1 or 27.k = 1 or 8, IB = 10to 18. For the top edges JB = NC-1 or 
NC+9 and KB = 1. For the bottom edges; one term will be JB = NC-1l or NG sain 
KB = 3, for the other term JB = NC or NC +8 and KB = 4. 

g. The Top and Bottom Side Edges (Medium Interface) 

This is the same as the end edges. The top equation has either the +1 or 
-1 index term replaced bv a medium grid difference term and the -k index term 1s replaced 
by the two surface terms for convection and radiation. The bottom equation has either 
the +1 or - j index term and the +k index term both replaced by a medium grid differ- 
ence term. These equations are valid over the range: 1 = 1 or 27,] = 2 to 26,k = 1 
or 8, JB = NC to NC +8. For the top edges IB = 9 or 19 and KB = |. For the bottom 
edges: one term will be IB = 9 or 19 and KB = 3, for the other term IB = 10 or 18 and 
KB = 4. 

h. The Corner Edges (Medium Interface) 

These have no surfaces nodes, just interfaces with the medium grid. The 

equations have both the +1 or -1 index term and the +j or -j index term replaced by a 


medium grid difference term. These equations are valid over the range: 1 = 1 or 27, ] 
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peor. kK — 2to 7, KB = 1to 3. One term will be IB = 9 or 19 and JB = NC or 
NC +8 and the other term will be 1B = 10 or 18 and JB = NC-1 or NC +9. 
1. Ihe Top and Bottom Corners (Medium Interface) 

Each of these eight equations is unique. They take into account a surface 
term for those on the top and an interface term for those on the bottom. In addition 
they each have two interface terms in the +1 or -i index and the +} or -} index. The fine 
grid nodal points that are corners are: (1,1,1) (1,1,8) (1,27,1) (1,27,8) (27,1,1) (27,1,8) 
ez) (27.27.8). 

3. Subroutine MED Fortran 

This subroutine contains the semi-discrete form of the explicit finite difference 
equation for the medium grid. The parameters for this grid are listed in Table 4. The 
subroutine 1s broken into thirteen sections, each dealing with a different tvpe of node. 
These nodes are shown in Figure 33 and Figure 34, where the letters correspond to the 
following sections. The constants in Table 4 are as normally defined for Fourier number 
and Biot number. Note that the Biot number Is with respect to Az, not Ax and that Az 
is one half of Ax for a surface node. 

The medium grid subroutine 1s the most complex because 1t must handle inter- 
facing with both the fine and the coarse grids. The distance to a coarse node 1s twice 
that of normal medium node spacing. Because the change in temperature at a node 1s 
related to 1/Ax, this reduces the effect by one half. With interfacing to the fine grid, the 
spacing is two-thirds of normal medium node spacing. In this case the effect is increased 
by three-halves. Interfacing with the fine grid also requires one additional step. Since 
one medium node is in contact with nine fine grid nodes (six for surface nodes), it is 
necessary to take the average temperature of the fine nodes. The second difficulty in 
interfacing is keeping track of the medium nodal points location with respect to the 
coarse and fine grids. In the subroutine all coarse array indexes are combined with an 
A to indicate they are A arrav subscripts and all fine array indexes are combined with a 
C to indicate they are C array subscripts. It 1s noted that the coarse, A, arrav 1s onlv 
two dimensional. 

a. The Interior Nodes 

The equation for the interior node is the simplest. This equation is valid 
over the following range: 1 = 2 to 26,j = 2 to 26 and k = 2 to 8 with the exception 
anv node adjacent to the fine zone, which must be handled separately. The zone formed 


by the nodes excluded is complicated and 1s as follows: fork < 4 theni = 10 to 18 and 


Table 4. PARAMETERS DEFINING THE MEDIUM GRID 


Dimensions: 81mm x 8Imm x 27mm 
Node Volume: 3mm x 3mm x 3mm 
Array: BU27710) 


The following constants are used: 
FO(2) = 1.4722 At 


B1(2) = 1.132x10° 
K = 53.0 W/mK 
Cp = 4.0x10° J/m?K 
h = 10.0 W)/m?K 





} = NC to NC +8, also fork < 4 with: = 9 or 19 and) = NC to NG+8, andialsamas, 
k <4withi= 10 to 18 andj = NC-1] or NC+9. 
b. The Top and Bottom Surface Nodes 

The surface nodes have the +k or -k index term replaced by a surface con- 
vection term. In addition, the top surface has an exclusion zone caused by the fine grid. 
The equations apply over the following range: 1 = 2 to 26,) = 2 to 26 and k = 1] or 
10. At the surface, k = 1, the exclusion zone 1s similar to the interior nodes: 1 = 10 to 
18 and} = NC to NC +8 also with 1 = 9 or 19 then}; = NC to NC +8, and with 1= 
1019 [8 then |= NN G-FomN G7! 

c. The Exterior End Faces (Coarse Interface) 

The end face equations have either a +1 or -i index term replaced with a 
difference term from the coarse zone. Due to the further distance of a coarse nodal point 
it is multiplied by one-half. These equations are valid over the following range: 1 = 2 
to 26,) = lor 27,k = 2to9, IA = 7to 15 and JA = NB-1 or NB+9. 

d. The Exterior Side Faces (Coarse Interface) 

The side face equations have either a +] or -} index term replaced with a 
difference term from the coarse zone. Due to the further distance of a coarse nodal point 
it is multiplied by one-half. These equations are valid over the following range: 1 = | 
or 27,}) = 2to 26,k = 2to9, IA = 6o0r l6and JA = NBto NBT8. 

e. The Exterior End Edges (Coarse Interface) 
The end edge equations have either a +1 or -1 index term replaced with a 


difference term from the coarse zone and a +k or -k index term replaced by a surface 
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Figure 33. The Medium Grid Node Types (Coarse Interface): The letters corre- 


spond the sections of the subroutine MED. 


convection term. These equations are valid over the following range: 1 = 2 to 26,j = 
momerok = | or 10, 1A = 7 to l5andJA = NB-! or NB+9. 
J. The Exterior Side Edges (Coarse Interface) 

The side edge equations have either a +j or -j index term replaced with a 
difference term froin the coarse zone and a +k or -k index term replaced by a surface 
convection term. These equations are valid over the following range: 1 = 1 or 27,j = 
2to 26,k = 1 or 10, IA = 6or 16andJA = NBto NB+S8. 

g. Ihe Exterior Corner Edges (Coarse Interface) 
The corner edges have both a +1 or -1 index and a +j or -j index terms re- 


placed with a difference term from the coarse zone. These equations are valid over the 
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Figure 34. The Medium Grid Node Type (Fine Interface): The figure is orientated 
upside-down to show the medium nodes which are adjacent to the fine 


zone. The letters correspond the sections of the subroutine MED. 


following range: 1 = 1 or 27,j] = 1 or 27,k = 2to9 and both 1) IA = 7 or I5 and JA 
= NB-I or NB+9 and 2) IA = 6 or 16 and JA = NB or NB+8. 
h. The Exterior Corners (Coarse Interface) 

The eight corners are all distinct equations. Thev each have two terms re- 
placed by a difference term from the coarse zone and one term replaced by a surface 
convection term. The equations are for the following nodal points: (1,1,1) (1,1,10) 
(1,270) (1,27, 10) (Q7 yee) (27527 ae 

1. Ihe Bottom of the Exclusion Zone (Fine Interface) 

For this equation the -k index term is replace with a difference term with 

nine fine zone nodes. This equation is valid for the following region: i = 10 to 18, j 


=NC to NCO =" =]"t38e 27 6 = eae 
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J. The Sides of the Exclusion Zone (Fine Interface) 

For this equation the +1 or -1 index term is replaced with a difference term 
with nine fine zone nodes. These equations are for the following region: 1 = 9 or 19, j 
= NC to NC+8,k = 2to3,IC = 1 or 27,JC = 1to KC = 3108. 

k. The Ends of the Exclusion Zone (Fine Interface) 

For this equation the +j or -j index term is replaced with a difference term 
with nine fine zone nodes. These equations are for the following region: 1 = 10 to 18, 
j = NC-1 or NC+9,k = 2103, IC = 1 to 27, JC = 1 or 27 and KC = 3to 8. 

l, The Side Edge of the Exclusion Zone (Fine Interface) 

For this equation the +1 or -i index term 1s replaced with a difference term 
with six fine nodes and a -k index term is replaced with a surface convection term. These 
equations are for the following 27 and region: 1 = 9 or 19,; = NC to NC+8,k = 1, 
IC = 1 or 27, JC = 1 to 27 and KC = | to 2. 

m. The End Edge of the Exclusion Zone (Fine Interface) 

For this equation the +] or -j index term is replaced with a difference term 
with six fine nodes and a -k index term 1s replaced with a surface convection term. These 
equations are for the following region: 1 = 10 to 18,} = NC-1 to NC+9,k = 1,1C 
= 110 27, JC = 1 or 27 and KC = | to 2. 

4. Subroutine COR Fortran 

This subroutine contains the semi-discrete form of the explicit finite difference 
equation for the coarse grid. The parameters for this grid are listed in Table 5. The 
subroutine is broken into six sections, each dealing with a different tvpe of node. These 
nodes are shown in Figure 35, where the letters correspond to the following sections. 
The constants in Table 5 are as normally defined for Fourier number and Biot number. 
Note the area of a node control volume surface is one-third the area of its side and that 
each node control volume has two surfaces. This modifies the Biot number to two-thirds 
of what it would be using a Ax spacing. Since the grid is two dimensional, all nodes have 
a surface convection term. 

The coarse grid subroutine is the simplest since it is two dimensional. On in- 
terfacing to the medium grid, the medium grid at the interface is two-thirds the distance 
of a normal coarse node. This increase the effect three-halves. Interfacing with the 
medium grid also requires one additional step. Since one coarse node is in contact with 
30 medium grid nodes it is necessary to take the average of the medium nodes. The 


second difficulty in interfacing is keeping track of the medium nodal points location with 


fo 


Table 5. PARAMETERS DEFINING THE COARSE GRID 


Dimensions: 180mm x 3l5mm x 27mm 
Node Volume: 9mm x 9mm x 2/mm 
Array: A(21,36) 
The following constants are used: 


FO(1) = 0.1635 Ar 


BI(1) = 1.132x107° 
K = 53.0 W/mK 
Cp = 4.0x10° Jim? K 
h = 10.0 W)/m?K 





respect to the coarse grid. In the subroutine all medium array indexes are combined with 
a B to indicate thev are B array subscripts. 
a. The Corners 
The selection criterion of the plate size was that the edges of the plate would 
not have a significant temperature change. Thus the four corners are assumed clamped 
and hence adiabatic. The indexes for the four corner are: (1,1) (1,36) (21,1) (21,36). 
b. The Ends 
It is assumed that negligible heat 1s lost out of the ends and hence either the 
+} or -j index term is dropped. These equations are valid over the following range: 1 
= 2to 20,} = 1 or 36. 
c. The Sides 
It is assumed that negligible heat 1s lost out of the sides and hence either the 
+1 or -} index term is dropped. These equations are valid over the following range: 1 
= Jor 2), = 210s. 
d. The Interior Nodes 
This is the basic equation and 1s valid over the following range: 1 = 2 to 
20, } = 2 to 35 with the exception of any node adjacent to the medium zone. The zone 
formed by the nodes excluded is as follows:1 = 7 to 15 and} = NB to NB+8, also with 
1= 6or 16 and} = NB to NBt+8, and also with i= 7 to 15 andj = NB-1 or NB+9. 
e. The Ends of the Exclusion Zone (Medium Interface) 
For this equation the +] or -] index term is replaced with a difference term 


with the 30 medium nodes. Care is taken to note that the surface medium nodes have 
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Figure 35. The Coarse Grid Node Types: The letters correspond to the subrou- 
fener © Te 


only have the volume of a regular medium node. These equations are valid over the 
range: 1 = 7 to 15,j) = NB-1i or NB+9, IB = 1 to 27, JB = 1 or 27 and KB = | to 
[0 
I. The Sides of the Exclusion Zone (Aledinin Interface) 
For this equation the +1 or -1 index term is replaced with a difference term 
with the 30 medium nodes. Care is taken to note that the surface medium nodes have 
onlv have the volume of a regular medium node. These equations are valid over the 


range: 1 = 6o0r 16,) = NB to NB+8, 1B = 1 or 27, JB = 1 to 2/7 and KB = 1 to 10. 


B. USER’S GUIDE TO MODIFYING AND RUNNING THE WELD PROGRAM 
The program WELD and its variants discussed in Appendix C are not interactive 


programs. Due to the amount of machine time that a typical simulation required, these 


8] 


programs were run on the grave shift by the batch processor, VMSCHED. Prior to each 
run, the program starting parameters were reviewed, and the program was recompiled. 
In Table 2 the major program variables are listed. There are several major changes that 
the user can decide to make before running the program: 

e The time for the program to run 

e Whether the run is new or a restart of a previous run 

e The arc velocity in the v and x directions 

e The arc power level 

e The starting temperatures 

e The material properties 

e The arc shape 

e Programming the arc movement and power history 


e Changing the grid sizes 


The first four of these are very simple to perform and are usuallv checked prior to each 
run. The last five are more complicated, and some guidance will be given about how to 
successfully change these parameters. 
1. Setting up for a normal run 
This section discusses how to set the time to run the program, the type of run, 
the arc power and speed. The time to run 1s specified by the variable FINI, just set this 
variable equal to the amount of time to simulate, in seconds. The next choice 1s to start 
a new simulation or restart an old one. To start a new simulation, insure that no old 
output files are on vour A disk, and comment out the statement “GOTO 100’. To restart 
a Old problem, insure that the output files are on your A disk and make the statement 
‘GOTO 100’ active. For problems which are not a restart you can also check the fol- 
lowing. The arc velocity is set bv the variable YVEL. This is the speed of the aroun 
millimeters per second. In the WELDMA version of the program, you can also set an 
x velocity, XVEL. However in the program WELD this variable is not used. The arc 
power level is the product of three variables; VOLT, AMP, and EFF. Where VOLT is 
in volts, AMP is in amperes, and EFF is a fraction of one, for the arc efficiency. This 
gives the power in watts. 
2. Setting the starting temperatures 
There are three places the starting temperatures are set. The first two are in the 


DATA statement. The A and B arrays (coarse and medium grids) are set by specifving 
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a temperature in degrees Kelvin. For the program listing below thev are set at 300 de- 
grees Kelvin. The starting temperature for the C array (fine grid) is specified bv its en- 
thalpy in watts per meter cubed. For the plain carbon steel being modeled, this is 
1.14237x10? HV/m3 at 300 degrees Kelvin. The last temperature that can be set is the 
temperature of the surroundings by the variable TINF. This also is in degree Kelvin. 

3. Changing the material properties 

In Table 3, Table 4, and Table 5, care was made to always define constants 
used for the Fourier and Biot numbers in each grid. In the program there is one block 
in Which these numbers are defined. By changing these constants, FO and BI, and set- 
ting the interface thermal conductivity, FKB, for the FIN subroutine, the material pro- 
perties and the surface conditions are changed. For the fine zone, which uses variable 
thermal properties, it is necessary to reprogram the three functions FK(T), T(H) and 
H(T). At present these are piece-wise continuous linear functions, but any function can 
be used. It is recommended that they be kept simple since they are repeatedly called by 
the program. 

When the material properties or mesh spacing are changed, the Fourier number 
is also changed. This mav require changing the step size to insure the stability require- 
ments are met, as discussed in Appendix B. When DELT is changed, it will also require 
changing the calculation of NDIV. These will be required to be changed not only at the 
Start of the program, but for restarts, in the restart section as well. If you will be 
changing DELT often, you may desire to change the program to calcuate NDIV as a 
function of DELT directly. 

4. Shaping the arc power distribution 

The program is written to allow shaping the arc as a double-elliptical gaussian 
distribution. This feature usually used a spherical gaussian distribution because arc 
shaping was not found to be necessary. Additionally it was feared that arc shaping 
would interfere with studving transient problems, since this would superimpose this em- 
pirical steady-state solution. The general form of the equation is given by equation A.2, 
Where the coefficients a, b, and c are the desired radial dimensions of the arc in the x, 
y, and z directions, respectively. The distances are in millimeters. In the source code, 
the use of different radii in the v direction has been commented out. This shows how 
the double elliptical pattern could be inputted for the arc power distribution. Note that 
it is required to change the source code twice. Once for calculating the volume averaging 
factor, SUM. And a second time for where the power is actually added to the nodes in 


the model. 
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5. Programming the arc movement and power history 

The arc movement is presently determined by two variables. YARC and DIS. 
Both of these variables are functions of YVEL* TIME. This product can be replace by 
anv desired function in both YARC and DIS defining equations. Thus arc accelerations, 
Starts and stops can be programmed in as desired. It is important to change both of 
these variables to keep track of the actual distance the arc has traveled. 

The arc power is in general constant and 1s inputted by the variable Q, which 
is the change of enthalpy of a standard node at one time step, that is in Joules per meter 
cubed. The easiest way to change Q is probably by multiplying it by a gain factor in the 
main DO loop. Thus, as the program runs you can change the arc power as desired. 
If the size of the grid spacing is being changed (i.e. from one millimeter spacing to 2 
millimeter spacing) then the variable Q would also be effected. For the example of two 
millimeters, then Q = QDOT*DELT*1.25E8. 

6. Changing the grid geometry 

To change the number of nodes is the most difficult change to perform. To 
change just the spacing of the grid, without changing the number of nodes, is just 
equivalent to changing the material properties. Just recalculate all of the material coef- 
fiecients, FO and BI. In addition vou will need to change the dimensions used in adding 
power; both for calculating the gaussian distribution as well as the change in enthalpy 
variable, Q. 

If the number and distribution of nodes is to be changed, than a major rewrite 
of the program 1s required. This would require checking every equation in the three su- 
broutines FIN, MED and COR as well as checking the main program sections which 
handle the input, output, Runge-Kutta sum, and grid shifting. An example of where this 
has been done is included in Appendix C. This 1s the version WELDC, where the fine 
and medium grids have been elongated to allow studing cooling rates. The modification 
approximately doubled the number of nodal points and hence doubled the run time of 


the program. 


C. SOURCE LISTING OF PROGRAM WELD FORTRAN 


PROGRAM WELD 


ts carhaccied ar beter ior bir birt oiek ep bak fae arr foro (ore oar oie pe oe Loe ck tok oars oy ok fo tok bok bk be io ied ar ak ord scr ier ok Lior dsr Orbe ier iwi ein oririscriseorisrisroar 
ve ve 
** DATE: 16 AUGUST 1988 WRITTEN BY: ROBERT ULE % 
ve ve 


a 2a a ae a ele a a a a a Fal a 
Yetedesesedetedeteseteteteletetiseieicicigigigiicnioicndkdilvvcddkvwcccediccceceki cock ichlesctek hicks Raced 


DIMENSION A(21,36),AOUT( 21,36) ,AIN(21, 36) ,ASUM( 21,36), 
*B(27,27,10),BOUT(27,27,10),BIN(27,27,10),BSUM(27,27,10), 
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e027 ,27,8),COUM 27,27,8),CIN(27,27,8),CSUM(27,27,8), 
*FO(3),BI(4) 

INTEGER NDIV,I,J,K,M,BSTEP,CSTEP 

CHARACTER ANS 

LOGICAL YES 

COMMON NB,NC,YARC,SUM,TINF,BI,FO 


DATA A,B/8046*300. 0/,C/5832"1. 14237E8/,ASUM, BSUM/8046*0. 0/ 
% AOUT , BOUT/8046*0. 0/ , CSUM/5832"0. / ,COUT/5832*0. / 


C SET TIME TO RUN THE PROBLEM AND FIXED CONDITIONS 


PINT = 5.00 
NDIV=FINI*100 
XVEL=. 0 
DELT=. 01 
TINF=300. 0 


C IF THIS IS A RESTART OF A PREVIOUS PROBLEM GOTO 100 
€ GOTO 100 
C SET INITIAL VALUES FOR STARTING A PROBLEM 
YVEL=4. 
VOLT=30. 
AMP=265. 
EFF=. 32 
C OPEN THE OUTPUT FILES 
OPEN(1,FILE='SURF' , STATUS='NEW' , FORM=' UNFORMATTED’ ) 
OPEN(2,FILE='FINAL ,STATUS='NEW ,FORM='UNFORMATTED' ) 
OPEN(3,FILE='CUT' ,STATUS='NEW' , FORM=' UNFORMATTED' ) 
C THE INITIAL CONDITIONS 
TIME=0. 
STEP=3. 
OUT=. 49 
N=0 
BSTEP=0 
CSTEP=0 
NB=3 
NC=10 
C WELD PARAMETERS 


QDOT=EFF*VOLT*AMP 
Q=QD0T*DELT*1. E9 


C REENTRY FOR RESTART 
200 CONTINUE 


C BOUNDARY CONDITION COEFFICIENTS 
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FOC 1)=DELT*. 1636 
FO(2)=DELT*1. 4722 
FO( 3 )=DELT*1000000. 
BIC i)=. 0022 
BI(2)=. 001132 
BI(3)=DELT*. 00009299 
BI(4)=DELT*50000. 


CALCULATE THE RUNGE-KUTTA APPROXIMATION 


DO 10 M=1,NDIV 
TIME=TIME+DELT 
DIS=TIME*YVEL 
N=N+1 


POSITION HEAT SOURCE AND CALCULATE VOLUME WEIGHTING FACTOR SUM’ 


YARC=YVEL*TIME+73-NB*9-NC*3 

LOF=65 - INT( VEL“ TIME) 

SUM=0. 

DO. 1s) 23 

IF ((J-YARC).GT.(0.0)) THEN “***** Allowed shaping the arc. 
SUM=SUM+QDENSE/4. *EXP( -. 10625*(( J-YARC )**2) ) 
ELSE 
SUM=SUM+QDENSE/10. *EXP( -. 017*(( J-YARC)**2) ) 
ENDIF 

lt CONTINUE 
SUM=SUM/Q 


ADD THE HEAT FROM THE ARC 


DO L2 =) e268 

IF ((J-YARC).GT.(€0.0)) THEN “***** Allowed shaping the arc 
Y=, 25*EXPC=. 10625( CJ-YARG)-2)) 

ELSE 
Y=. 1“EXP( -. 017*C(C J-YARC)**2) ) 

ENDIF 

DO 2 I=11,17 
DO 2 kK=1,3 

CCl,J,K )=CC 1.0, K )4AY*32(¢ 1-105) Sun 
2 CONTINUE 


CONVERT ENTHALPY TO TEMPERATURE FOR THE FINE PLATE 


DO 60 I=1,27 
DO 60 J=1,27 
‘DO 60 K=1,8 
CIN(I,J,K)=T(C(1I,J,K)) 
60 CONTINUE 


FIRST RUNGE-KUTTA ITERATION 
CALL COR(A,B,ASUM) 


CALL MED(A,B,CIN,BSUM) 
CALL FIN(B,CIN,CSUM) 
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FIND T+K1/2 


oat I=. 01 
DOr i= i0e6 
11 AIN(I,J)=A(1,J)+ASUM(I,J)/2. 
DO 12 I=1,27 
DO 12 J=1,27 
DO 13 K=1,10 
13. BIN(I,J,K)=B(1,J,K)+BSUM(1I,J,K)/2. 
Dome K-1,8 
CIN(1,J,K)=T(C(1,J,K)+CSUM(1I,J,K)/2.) 
12 CONTINUE 


SECOND ITERATION 


CALL COR(AIN,BIN,AOUT) 
CALL MED(AIN,BIN,CIN, BOUT) 
CALL FIN(BIN,CIN, COUT) 


FIND T+K2/2 


DO 21 I=1,21 
DO 21 J=1,36 
ASUM(I,J)=ASUM(I,J)+2. *AOUT(I,J) 
21 AIN(I,J)=A(I,J)+AOUT(I,J)/2. 
Den?) 1=1,27 
DO 22 J=1,27 
DO 23 K=1,10 
BSUM(I,J,K)=BSUM(I,J,K)+2. *BOUT(I,J,K) 
23. BIN(I,J,K)=B(1I,J,K)+BOUT(1I,J,K)/2. 
DO 22 K=1,8 
CSUM(I,J,K)=CSUM(I,J,K)+2. *COUT(I,J,K) 
CIN(I,J,K)=T(C(1,J,K)+COUT(I,J,K)/2. ) 
22 CONTINUE 


THIRD ITERATION 
CALL COR(AIN,BIN,AOUT) 


CALL MED(AIN,BIN,CIN, BOUT) 
CALL FIN(CBIN,CIN,COUT) 


END f+K3 
DO 31 [=1,21 
Oe sted -1, 36 


ASUM(I,J)=ASUM(I,J)+2. *AOUT(I,J) 
31 AIN(I,J)=A(I,J)+AOUT(I,J) 
Dons? I=1,27 
DONS? J=1,27 
DO) 33 K=1/10 
BSUM(I,J,K)=BSUM(I,J,K)+2. *BOUT(1I,J,K) 
33. BIN(I,J,K)=B(I,J,K)+BOUT(I,J,K) 
DO 32 K=1,8 
CSUM(I,J,K)=CSUM(1I,J,K)+2. *COUT(I,J,K) 
CIN(I,J,K)=T(C(I,J,K)+COUT(I,J,K)) 
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32 CONTINUE 
FORTH ITERATION 


CALL COR(AIN, BIN, AOUT) 
CALL MED(AIN, BIN,CIN, BOUT) 
CALL FIN(BIN,CIN,COUT) 


PERFORM THE RUNGE-KUTTA SUM 


0) (AG, ea 21 
DO 40 J=1,36 
40 ACI,J)=A(1I,J)+( ASUM(I,J)+AOUT(I,J))/6 
DO 41 I=1,27 
DO 41 J=1,27 
DO 43 K=1,10 
43°  B(I,J,K)=B(1,J,K)+(BSUM(I,J,K)+BOUT(I,J,K))/6 
DO 41 K=1,8 
41 €(€1;J,K)=CC1 J. Ky4CCSUM(1, J, K)+COUT(C TK) 76 


STEPSGRID IF sREQUTRED 


119 (OES Ea, SIMA) Galen 
STEP=STEP+3. 
BSTEP=BSTEP+1 


MOVE FINE GRID, FIRST AVERAGE FINE TEMPS AND ASSIGN TO MED GRID 


DO 51 IB=1,9 
TB1=0. 
TB2=0. 
TB3=0. 
DO 50 IC=IB*3-2,IB*3 
DOr sSeamiG= 155 
TRI=T(COIC JC, 1) 42. CeCe JG72) 4 nbd 
DO 50 KC=3,5 
TB2=TB2+T(C(IC,JC,KC)) 
TB3=TB3+T(C(IC,JC,KC+3)) 
50 CONTINUE 
BC IB+9,NC,1)=TB1/27 
BC IB+9 ,NC,2)=TB2/27 
B( IB+9 ,NC,3)=TB3/27 
51 CONTINUE 


SHIFT FINE GRID 
BO) 32 Je we 
JJ=J+3 
DO 52 I=1,27 
DO 52 K=1,8 
52 CCl ak) =CC cree) 
ASSIGN MED TEMPS TO FINE GRID 


DO 53 I=10,18 
TB1=BCI,NC+9,1) 
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TB2=B( I ,NC+9,2) 
TB3=B(I,NC+9,3) 
DO 53 IC=(1-9)*3-2,(1-9)*3 
DO 53 JC=25,27 
CC1IG, IC .1)=H(TB1) 
GONG. 0G. 2) =—HCne 1) 
DO 53 KC=3,5 
CCUG, WOO Sa 
C(IC,JC ,KC+3)=H(TB3) 
53. CONTINUE 
NC=NC+1 


Cinch 10 SER IF TIME TO SHIFT MEDIUM GRID 
PP UBeatEP. BQ. 3) THEN 
BSTEP=0 


THEN AVERAGE MEDIUM GRID AND ASSIGN TO COARSE GRID 


DO 54 IA=1,9 
TA1=0. 
DO 55 IB=IA*3-2,IA*3 
DO 55 JB=1,3 
TA1=TA1+( B(IB,JB,1)+B(IB,JB,10))/2 
DO 55 KB=2,9 
TA1=TA1+B(IB,JB,KB) 
55 CONTINUE 
A( IA+6 ,NB)=TA1/81 
54 CONTINUE 


SHIFT MED GRID 


DO 56 J=1,24 
JJ=HI+3 
DO 56 I=1,27 
DO 56 K=1,10 
56 B@led 1) Bl sda) 


ASSIGN COARSE TEMPS TO MEDIUM GRID 


DO 57 I=7,15 
TA1=A(I,NB+9) 
DO 57 IB=(I-6)*3-2,(1-6)*3 
DO 57 JB=25,27 
DO 57 KB=1,10 
B(IB,JB,KB)=TA1 
57 CONTINUE 
NB=NB+1 
NC=10 
ENDIF 
ENDIF 


OUTPUT SOLUTION FILES 


Pi ClLIME,GE.OUT) THEN 
OUT=OUT+. 5 
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WRITE( 1) (CTUCC(1 ai), 1=127 ) =e 
C AVERAGE FINE TEMPS AND ASSIGN TO MEDIUM GRID 
DO 61 J=0,8 
DO 61 IB=1,9 
TB1=0. 
DO 66 IC=IB*3-2,1B*3 
DO 66 JC=J*3+1,JI*3+3 
66 TBI=TCEC Te; Je,1)) +271 CCl, IC.2 eames 
61 B(9+1B,NC+J,1)=TB1/27. 
WRITE(1)((BC1,J,1),1=1,27),J=1,27) 
WRITE( 3)((TC(C(I, INT(YARC) ,K)),I=1,27) ,K=1,8) 
ENDIF 
IF (NB. LE. 11) THEN 
IB=26-(NB-3)*3 
IF (IB. GE. NC. AND. IB. LE.NC+8) THEN 
IC=2+( IB-NC)**3 
WRITE(4) (T(C(IC,J,1)),J=2,26,3) 
ELSE 
WRITE(4) (B(IB,J,1),J=10,18) 
ENDIF 
ENDIF 


10 CONTINUE 
C OUTPUT FINAL RESULTS 


WRITE(2) TIME 
WRITE(2)(((T(C(1,J,K)),I=1,27),J=1,27) ,K=1,8) 
C AVERAGE FINE TEMPERATURES AND ASSIGN TO MEDIUM GRID 
DO 62 J=0,8 
DO 62 IB=1,9 
TB1=0. 
TB2=0. 
TB3=0. 
DO 63 IC=IB*3-2,1B*3 
DO 63 JC=1+J5*3 ,3+I*3 
TBI=T(C(1IC,JC, 1))+2087(CC(IC,IC,2))4+7 BL 
DO 63 KC=3,5 
Reo —h 2 teu lee meen 
TB3=TB3+T(C( IC, JC,KC+3)) 
63 CONTINUE 
B(IB+9 ,NC+J, 1)=TB1/27. 
BC IB+9 ,NC+J , 2)=TB2/27. 
B( IB+9 ,NC+J ,3)=TB3/27. 
62 CONTINUE 
WRITE(2)(((B(1,J,K), 1=1,27)),J=1,27) ,k=aelon 
C AVERAGE MEDIUM TEMPERATURES AND ASSIGN TO COARSE GRID 
DO 64 J=0,8 
DO 64 IA=1,9 
TA1=0. 
DO 65 IB=IA*3-2,IA*3 
DO 65 JB=1+J*3,3+J*3 
TA1=TA1+(B(IB,JB,1)+B(1IB,JB,10))/2 
DO 65 KB=2,9 
TA1=TA1+B(IB,JB,KB) 
65 CONTINUE 
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AC IA+6 ,NB+J)=TA1/81 
64 CONTINUE 
Mol EEC ZC (ACI J), 1=1,21) ,J=1, 36) 


C ALLOW OPTION OF RESTARTING THE PROBLEM 


INQUIRE( 9 , OPENED=YES ) 

IF(YES) THEN 

REWIND(9) 

ELSE 

OPEN(9 ,FILE='RESTAR' ,STATUS='NEW' , FORM='UNFORMATTED’ ) 
ENDIF 
WRITE(9)OUT,Q,NB,NC,YVEL,STEP,CSTEP,BSTEP,N,DELT 
WRITE(9)A 
WRITE(9)B 
WRITE(9)C 

STOP 


C PROCEDURE FOR RESTARTING A PREVIOUS PROBLEM 
100 OPEN(2,FILE=' FINAL’ ,STATUS=' OLD’ , FORM='UNFORMATTED' } 
C THE INITIAL CONDITIONS 


READ( 2)}TIME 

REWIND(2) 

IF (TIME. LT. (.49)) THEN 
OPEN(1,FILE='SURF’ ,STATUS='NEW' , FORM='UNFORMATTED' ) 
OPEN(3,FILE='CUT' , STATUS='NEW’ , FORM='UNFORMATTED’ } 

ELSE 
OPEN(1,FILE=' SURF’ , STATUS='OLD' , FORM=' UNFORMATTED’ ) 
OPEN(3,FILE='CUT’ ,STATUS='OLD' , FORM='UNFORMATTED' } 

ENDIF 

OPEN(9,FILE='RESTAR’ ,STATUS=' OLD’ , FORM=' UNFORMATTED’ ) 

READ( 9 )OUT,Q,NB,NC,YVEL,STEP,CSTEP,BSTEP ,N,DELT 


feewoe LON THE RUNNING OUTPUT FILES TO THE END OF THE FILES 


DO 102 L=1,N/50 
READCIDCCECI (J,1)51 
READCIYCCBCI,J,1), 1 
ReApCSy (CCCI, !1,K),1 

102 CONTINUE 


Hou il 
+ b+ b+ 
NO BRO PO 
NINN 
Ne ee 
AGG 
Wout ll 
+ b+ 
OO RO ho 
ee 
w/ N/ 


C INPUT THE TEMPERATURE/ENTHALPY MATRICES 


READ(9)A 
READ(9)B 
READ(9)C 
GOTO 200 
END 


C FUNCTION FOR K AS A VARIABLE OF TEMPERATURE 


REAL FUNCTION FK(T) 
ell. 975, 2 TREN 


: 9] 


FK=63. 92-. 04*T 
ELSE 
TE Ciel 73.) HN 
PK=165 908/54, 00625-1 


C FUNCTION FOR FINDING TEMPERATURE AS A FUNCTION OF ENTHALPY 


REAL FUNCTION TCH) 
IP CHAE 2a or?) THEN 
T2277 3636b=7 a2 3 
ELSE 
EF (Ho bis 2259 EN 
MS 10 Ss ole a (215) 
ELSE 
TF CH. Diy SE eenHEN 
P1625 56-7 tooo 
ELSE 
TEC Or 5E9 ) THEN 
T=1. 8182E-8*H+1586. 64 
ELSE 
T=2. 222E-7*H-504. 45 
ENDIF 
ENDIF 
ENDIF 
ENDIF 
RETURN 
END 


GC FUNCTION FOR FINDING ENTHALPY AS A FUNCTION OF TEMPERATURE 


REAL FUNCTION H(T) 
IF Ch i 9225) THEN 
H=4. 231E6**T-1. 15506E9 
ELSE 
IF (T. LT. 948.) THEN 
H=. 02E9*T-15. 71E9 
ELSE 
IF (T. LT. 1723.) THEN 
H=5. 484E6*"T-1. 9488E9 
ELSE 
IF (T. LT. 1773.) THEN 
H=. 055E9*T-87. 265E9 
ELSE 
H=. 0045E9*T+2. 2715E9 
ENDIF 
ENDIF 
ENDIF 
ENDIF 
END 
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D. SOURCE LISTING OF SUBROUTINE FIN FORTRAN 
SUBROUTINE FIN(B,C,COUT) 


DIMENSUON C2727 10) C(27, 27 supe count 2742725) ,FOC3) BIC) 
SOUMON NB YNG, LARC, SUM, TINE ,BIGFO 


GK(T1,T2)=FK(T1)*(T1-T2) /(PR(T1)4FK(T2)) 
DATA FKB/26. 5/ 
* A. THE INTERIOR NODES 


Mon 1 1=2,26 
DO 1 J=2,26 
DOr ts Kaa 
COUTCI. JR )=hK(C( Je) OCS GRC C(14+1,J,K),CC1,J.K))+ 
Che Clete ek Clee) rc CG les ak cure) 4CKCCGl sJ-1.K), 
eC 1 JK) )4GR(C( 1.3, kt1) ,6(1,3,8) 4GKCCCl, J,K-1).CC1,J,K))) 
1 CONTINUE 


feo. tHE TOP SURFACE NODES 


TINF4=TINF"4 

DO 2 1=2,26 

DO 2 J=2,26 

SCCReIN 1) —MMCCCl 3,1) hOtCsn-CGRCCCi+ i J, 1),CC1,J,1))4+ 

PCKGe@re ls J-1) COL J, 1))+GRCGCI. gt) 1) CCl. 7, 1+ 
MCKGCt = tb) OC.) 2. “OKCCCl dee), CCl. 51) js 

*  BI(3)*(TINF4-C(1,J,1)**4)+B1I(4)**(TINF-C(1I,J,1)) 
2 CONTINUE 


* C. THE BOTTOM SURFACE NODES (MEDIUM INTERFACE) 


DO 3 I=2,26 
IB=(I-.5)/3+10 
DO 3 J=2,26 
JB=(J-.5)/3+NC 
GOUT(1,J,8)=FO(3)*(FK(C(1,J,8))*(GK(C(I+1,J,8),C(1,J,8))+ 
* GK(C(I-1,J,8),C(1,J,8))+GK(C(I,J+1,8),C(1,J,8))+GK(C(I,J-1,8), 
PC 1d 8) aC CCClss 7 (ee) Fhe BC 18. J8,4)-C(1.J,8))) 
3 CONTINUE 


* D. THE END FACES (MEDIUM INTERFACE) 


JB=NC-1 
JMAX=NC+9 
DO 4 I=2,27 
IB=(I-.5)/3+10 
DO 4 K=2,7 
KB=2 
IF(K.EQ.2) KB=1 
IF(K.GE.6) KB=3 
COUT(I ,1,K)=FO(3)*(FK(C(1,1,K))*(GK(C(1+1,1,K),C(1,1,K))+ 
CRC CGE ar GCl 1K) 4GK(C( Ll. 20K GCI wile) GKCCCL.1,K+1), 
*  C(1I,1,K))+GK(C(1,1,K-1),C(1,1,K)))+FKB*(B(1IB,JB,KB)- 
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GU) 3) 
COUT(I ,27 ,K)=FO(3)*(FKC(CC1,27,K) )*(GK(C(1+1,27,K),C( 1,27 seeies 
* GK(C(1I-1,27¢K)eGC1 ,27 ,K) )+GK(C( 1, Zope Ge kee 
*  GK(C(1,27,K4t1) 001, 27,K) )tGK(C( 1, 27 hese 1 9 atone 
* = FKB*(B( IB, JMAX,KB)-C(1I,27,K))) 
4 CONTINUE 


* E. THE SIDE FACES (MEDIUM INTERFACE) 


DO 5 J=2,26 
JB=(I-.5)/3+NC 
DOS k= 27 
KB=2 
ihCKeHOR2)ee— 
IF(K.GE.6) KB=3 
COUT( 1,J,K)=FO(3)*(FK(C(1,J,K))*(GK(C(2.J.K).C@ied shone 
*  GK(C(1, Jt1.K),6( 1,04) )+-GRCCC1 I= Knee eke) 
*  GK(C(1,J,K+1),C(1,J,K))+GK(C(1,J,K-1),C(1,J,K)))+ 
*  FKB*(B(9,JB,KB)-C(1,J,K))) 
COUT( 27 ,J,K)=FO(3)*(FK(C(27,J,K))*(GK(C(26,J,K) ,C(27,J,K))+ 
* GK(C(27,J+1,K),C(27,J,K) )+GK(C(27,J-1,K),C(27,J,K))+ 
*  GK(C(27,J,K+1),6(27,J,K) )+GK(C( 27, J. kat) 6 Ger kone 
*  FKB*(B(19,JB,KB)-C(27,J,K))) 
5 CONTINUE 


* F. THE TOP AND BOTTOM END EDGES (MEDIUM INTERFACE) 


KB=1 
JB=NC-1 
JMAX=NC+9 
DOG, I=2 26 
iB=(1-. 5) /sr10 
COUT(I,1,1)=FOC3)*CEKGCC1,1,1))* (GRC Iti ea 
: GK(CCI<1 tet Tere 2 earn eae 
* 2. #GK(G(T,1,2) CCI, 1,1)))4FRB<(BGLB, JB RB EOGmMNip phe 
* BIC 3)*CTINF4-C(1I,1,1)**4)+B1(¢4)*( TINF=-C(1,1,1)) 
COUT(1,27,1)=FO(3)*(FK(C(1,27,1))*(GK(C(I+1,27,1),C(1,27,1))+ 
* GKC CC l-1,27 51),CC1,27,1))+GKCCCL. 265.1) Cela ala 
* 2.*GK(C(1,27,2),CC1,27,1)) )+EKB*( BCIB ,JMAX KE) =CC I jens 
* BIC 3)*(TINF4-C(1,27,1)**4)+BI(4)*(TINF-C(I,27,1)) 
COUTC1I,1,8)=FO(3)*(FKCCC1,1,8))*(GKCCC 11 3) CCl, 1, Soe 
os GRhCC Olen 1, 1, 5) FC NCC 2 Greeley scr 
*  GK(C(I,1,7),C(I,1,8)))+FKB*(B(1B,JB,3)+B( IB,NC,4) 
aS pre OlG Ieee 169) 
COUT(1,27,8)=FO(3)*( FK(C(1,27,8))*(GK(C(I+1,27,8),C(I,27,8))+ 
*  GK(C(I-1,27,8),C(L,27,8) )+GK(C(1,26,8),C(1,27,8))+ 
*  GK(C(I,27,7),C(1I,27,8)))+FKB*(B( 1B, JMAX, 3)+B(IB,NC+8,4)- 
* 2.*C(1,27,8))) 
6 CONTINUE 


* G. THE TOP AND BOTTOM SIDE EDGES (MEDIUM INTERFACE) 
DO 7 J=2,26 
JB=CI-. 5) 755NC 


COUT(1,J, 1)=FOC3)*(FKCCC1, 3,1) )=(GRCCOIe I-41 1) CCl Is 
*  GK(C(1,J-1, 1), CCWs 1) )4CKhCC (2 Jd Cle ieee 
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* 2. GRC. da 2G) DEE BBC oes.) -C( lem. 1) ))s 
* BI(3)*(TINF4-C(1,J,1)**4)+BI(4)*(TINF-C(1,J,1)) 
COUTGZ TD) Fh OCS)CRRCECZ7 - al es(GhCCC2 7 Jl, 1), CC 27 ,J51))+ 
pee GN(O(27 J) CC 27 1) ) eGKeC@24 7 e027. 3.1))4 
% Pecan G@2 7.0.2) ,.C( 27, J. 1) )eEKB*( BC 1998 1) =C(27,J,1)))+ 
* BI(3)*(TINF4-C(27,J,1)**4)+BI(4)*(TINF-C(27,J,1)) 
COUT(1,J,8)=FO(3)*(FK(C(1,J,8))*(GK(C(1,J+1,8),C(1,J,8))+ 
oy CRUG JigleoeCCladne  reKC COZ ene la) .6))) + 
* GK(C(1,J,7),€(1,J5,8) ))+FKB*(BC9 ,JB,3)+B(10,JB,4) 
* § =2.*C(1,3,8))) 
COUT(27,J,8)=FO(3)*( FK(C(27,J,8))*(GK(C(27,J+1,8) ,C(27,J,8))+ 
* GK(€(27,J-1,8) ,C(27,J,8))+GK(C(26,3,8),€(27,J3,8))+ 
* GK(C(27,J3,7),C(27,J,8) ) )+FKB*(B(19,JB,3)+B(18,JB,4) 
¥ ee C2 7G) i) 
7 CONTINUE 


* H. THE CORNER EDGES (MEDIUM INTERFACE) 


DO 8 K=2,7 
KB=2 
IF(K. EQ. 2) KB=1 
IF(K.GE.6) KB=3 
Oba @ineieG)—FOGs)GhK(C Gilet. K peeGR(Gr 2, 1.K).c(1,1,K))+ 
oC c kh) Ol, Ko) rGKCe, 1) Thre Gla kn 
*  GK(C(1,1,K-1),C(1,1,K)))+FKB*(B(10,NC-1,KB)+B(9,NC,KB)- 
me 22°06 1.1,K))) 
mOUM@ne 7 kh \=BOCS a (PKCCCL, 27K) (GRCC(2.27.K).C(1,.27.K)4 
PGE GCIs 26nk ) OClN27 ,.K) 4GKCC Cine 7 Kd eG 27K) 
*  GK(C(1,27,K-1),C(1,27,K)) )+FKB*(B(10,NC+9 ,KB)+B(9,NC+8,KB)- 
sae? 01 ,27,K))) 
GOUT( 27, 1,K)=FO( 3)**( FK(C( 27, 1,K))*(GK(C( 26, 1,K).C(27,1,K))+ 
Ga GRC 279 -K)..60 271. KY) +GKCCC27 41 Kt) -CO27e1.K))+ 
*  GK(C(27,1,K-1),C(27,1,K)))+FKB*(B(18,NC-1,KB)+B(19,NC,KB)- 
nae 2 °0(27 ,1,K))) 
OUM@2e 2 7K )=FOCs)*( FKCCC27 27 .K) )=(GKCG( 27, 27 eK) VC027 .27,.K))+ 
mae GER 6( 27, 26.K),C( 27, 27.K))4+GK(C( 27.27 ,.K+1),C(27,27,K))+ 
* GK(C(27,27,K-1) ,C(27,27,K)) )+FKB*(B(19,NC+9 ,KB)+ 
*% BC 18,NC+8,KB)-2.*C(27,27,K))) 
8 CONTINUE 


* I. THE TOP AND BOTTOM CORNERS (MEDIUM INTERFACE) 


eOURGl 1, 1)=FOC 3)*(PKCCC1. 1. 1) )*CGR(C0 2.1, 1)601,1,1))+ 
re COGin, 2esl CC 1 1) HGR CGI ie) CCl «1 tEKBy 
*  (B(10,NC-1,1)+B(9,NC,1)-2.*C(1,1,1)))+ 
% BI(3)*(TINF4-C(1,1,1)%*4)+B1(4)%(TINF-C(1,1,1)) 
BOURCIN7 1) =FOC3)7( PROOCI. 275 1 a (GK( CC. 27, 1).6(1,27,1))+ 
PGR Cal 2Gn 1) CC le 27) GRC OC lego) CC 127 1) ANKE 
PG ONG 10 14809 NC18e1)-22"c( 1.27.1) ))+ 
% 8I(3)*(TINF4-C(1,27,1)**4)+B1(4)*(TINF-C(1,27,1)) 
BOUme yam )—FOCl rt ( FK( C027, 1. 1 ))-(GKCCC26. 141) ,0(27,151))+ 
CO. 7, 2), CC 27 GR CG( 27, 1, 2), CC 27s) ERB 
ee Gbalemnc= 1 1) est 19.NC,1)-2. 0027, 1.1) )4+ 
* BI(3)*(TINF4-C(27,1,1)**4)+B1(4)*(TINF-C(27,1,1)) 
BOUMG2 78) =FOCs (ERCCC27.27 11) )=(GK0C(26.27,1),0027,27,1))+ 
Pea 7 2601) C027 27. 1) CRC 27.272) 027,27, 1)) )tEKB™ 
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* (B(18,NC+9 ,1)+B(19,NC+8,1)-2. *C(27,27,1)))+ 
* = BI(3)*(TINF4-C(27,27,1)**4)+BI(4)*( TINF-C( 27 ,27,1)) 
COUT(1,1,8)=FO(3)*(FK(C(1,1,8))*(GKCC(2,1,8).6(1,1,8))4 
* GK(C(1,2,8),CCla iG) CK CCI 7 eect te 
* (B(10,NC-1,3)+B(9,NC,3)+B(10,NC,4)-3.*C(1,1,8))) 
COUT( 1,27,8)=FO(3)*(FK(C(1,27,8) )*(GK(G( 2.27.8) G01, 27menier 
*%  GK(C(1,26,8).,G(1.27.8) )+GK(G(1,27,7).0(1,27.6)) ) +e 
* (BC 10,NC+9,3)+B(9,NC+8, 3)+B(10,NC+8,4)-3. *C(1,27,8))) 
COUT( 27, 1,8)=FO( 3)*(FK(C( 27, 1,8) )*(GK(C( 26,1,8) ,C(27.1,8))+ 
* GK(C(27,2,8),C(27,1,8))+GK(C(27,1,7),C(27,1,8)) )+FKB* 
*  (B(18,NC-1,3)+B(19,NC,3)+B(18,NC,4)-3. *C(27,1,8))) 
COUT( 27 ,27,8)=FO( 3)*(FK(C(27,27,8) )*(GK(C( 26, 27,8) ,C(27,27,8))+ 
*  GK(C(27, 26,8) ,C(27,27,8) )AGK(G(27 27,7 InGl 27.278 eee 
%  (B(18,NC+9,3)+B(19,NC+8,3)+B(18,NC+8,4)-3. *C(27,27,8))) 


RETURN 
END 


E. SOURCE LISTING OF SUBROUTINE MED FORTRAN 


SUBROUTINE MED(A,B,C, BOUT) 


DIMENSION A(21,36),BOUT(27,27,10),B(27,27,10),BSUM(27,27,10), 
*C(27,27,8),FO(3),B1(4) 
COMMON NB,NC,YARC,SUM,TINF,BI,FO 


*« <A. THE SINTERTOR NODES 


DO 7 I=2,26 
DO 7 J=2,26 
DO 777 K=2,4 
IF((I.GE. 9. AND. I. LE. 19. AND. J. GE. NC-1. AND. J. LE. NC+9) 
% | AND. . NOT. (K. EQ. 4. AND. (I. EQ. 9. OR. I. EQ. 19. OR. J. EQ. NC-1. OR. J. EQ. 
NC+9)). AND. . NOT. ((J. EQ. NC-1. OR. J. EQ. NC+9). AND. (I. EQ. 9. OR. I. EQ. 
1) 1G mGOnO 777 
BOUT(I,J,K)=FO(2)*( B( I+1,J,K)+B(1-1,J,K)+B(1I,J-1,K)+B(1,J+1,K)+ 
ve B(1.J,K+1)+B( 1.) ,K-1)-6) “PC. > 
777 CONTINUE 
DO 7 K=5,9 
BOUT(I,J,K)=FO(2)*(B(I+1,J,K)+B(1I-1,J,K)+B(1,J-1,K)+B(1,J+1,K)+ 
* B(I,J,K+1)+B(1,J,K-1)-6. *B(I,J,K)) 
7 CONTINUE 


* B. THE TOP AND BOTTOM SURFACE NODES 


DO 6 I=2,26 
DO 8 J=2,26 
BOUT(1I,J,10)=FO(2)*(B( I+1,J,10)+B(I-1,J,10)+B(1,J+1, 10)+ 
B(I,J-1,10)+2*B(1,J,9)+BI(2)*TINF-(BI(2)+6. )*B(1I,J,10)) 


IF((1. GE. 9. AND. I. LE. 19. AND. J. GE. NC-1. AND. J. LE. NC+9) 
* AND. .NOT. ((J. EQ. NC-1. OR. J. EQ. NC+9). AND. (I. EQ. 9. OR. I. EQ. 
* 19))) GOTO 8 
BOUT(1,J,1)=FO(2)*(B( I+1,J,1)+B(I-1,J,1)+B(1,J+1,1)+B(I,J-1,1) 
* +2*B(1I,J,2)+BI(2)*TINF-(BI(2)+6. )*B(I,J,1)) 
8 CONTINUE 
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feo. THE EXTERIOR END FACES (COARSE INTERFACE ) 


DO 9 I=2,26 
IA=INT((I-.5)/3)+7 
DO 9 K=2,9 
BOUT(I,1,K)=FO( 2)*(B(I+1,1,K)+B(I-1,1,K)+B(1,2,K)+B(I,1,K+1)+ 
B(I,1,K-1)+. 5*A( IA,NB-1)-5. 5°**B(I,1,K)) 
BOUT(1,27,K)=FO(2)*(BCI+1,27,K)+B(1-1,27,K)+B(1,26,K)+ 
ve B(1,27,K+1)+B(1,27,K-1)+. 5*A( IA, NB+9)-5. 5*B(1,27,K)) 
9 CONTINUE 


* D. THE EXTERIOR SIDE FACES (COARSE INTERFACE) 


DO 10 J=2,26 
JA=INT((J-.5)/3)+NB 
io 100 kK=2 9 
BOUT(1,J,K)=FO(2)*(B(2,J,K)+B(1,J+1,K)+B(1,J-1,K)+B(1,J,K+1)+ 
% B(1,J,K-1)+. 5%*A(6,JA)-5. 5*B(1,J,K)) 
BOUT( 27,J,K)=FO(2)*(B(26,J,K)+B(27,J+1,K)+B(27,J-1,K)+ 
B(27,J,K+1)+B(27,J,K-1)+. 5*A(16,JA)-5. 5*B(27,J,K)) 
10 CONTINUE 


* E. THE EXTERIOR END EDGES (COARSE INTERFACE) 
WO 11 I=2,26 


IA=INT((CI-. 5)/3)+7 
BOULCIY1, 1)=FO(2)* AC Glick ala si CLS alesse Gli Amo Ee a2: Ga lena ee 


ze 5%A( 1A,NB-1)+BI(2)#TINF-(4. 5+B1(2))*B(1I,1,1)) 
BOUT(I,1, 10)=FO(2)* oG@E eine IO) BGI iy lO) an, 10)+B(1,1,9)+ 
zs 5A(1A,NB-1)+B1(2)*TINF-(4. 5+BI(2))*B(1,1,10)) 
BOUT(I,27,1)=FO(2)**(B(I+1,27,1)+B(1-1,27,1)+B(1,26,1)+B(1,27,2)+ 
ve SAC IA,NB+9)+B1(2)"TINF-(4. 5+B1(2))*B(1,27,1)) 
BOUT(I ,27,10)=FO(2)*(B(I+1,27,10)+B(1I-1,27,10)+B(1,26,10)+ 
zs B(I,27,9)+. 5*A( IA,NB+9 )+B1(2)*TINF-(4. 5+BI(2))*B(1,27,10)) 


11 CONTINUE 
* F, THE COARSE SIDE EDGES (COARSE INTERFACE) 


BO 12 J=2,26 
JA=INT((J-. 5)/3)+NB 
BOUT(1,J,1)=FO(2)*(B(2,J,1)+B(1,J+1,1)+B(1,J-1,1)+B(1,J,2)+ 
* ,5*A(6,JA)+B1(2)*TINF-(4. 5+BI(2))*B(1,J,1)) 
BOUT( 27 ,J,1)=FO(2)*(B( 26 ,J,1)+B(27,J+1,1)+B(27,J-1,1)+B(27,J,2)+ 
ce -5%A(16,JA)+B1(2)*TINF-(4. 5+BI(2))*B(27,J,1)) 
BOUT(1,J,10)=FO(2)*(B(2,J,10)+B(1,J+1,10)+B(1,J-1,10)+B(1,J,9)+ 
ve .5%A(6,JA)+B1(2)*TINF-(4. 5+B1(2))*B(1,J,10)) 
BOUT(27,J,10)=FO(2)**(B(26,J,10)+B(27,J+1,10)+B(27,J-1,10)+ 
ve B(27,J,9)+. 5*A( 16, JA)+B1(2)*TINF-(4. 5+BI(2))*B(27,J,10)) 
12 CONTINUE 


* G. THE EXTERIOR CORNER EDGES (COARSE INTERFACE) 
HON 13 K=29 
, BOUT 1, ley =FOC2)*CB( 2, 05K) re lz »K)+BC1, Ke) Bl eeeen 


mac CACO Nema? yNB=1))iroe ee ld -K)) 
BOUTCZ7, | K)=FO(2)? eB CZGn re G27 2 Rte 1 Kd 4B027 1, K-1)4+ 


at 


¥ . 9*(AC16,NBICAGIS INB-1))=Sa iy K)) 
BOUT( 1,27 ,K)=FO(2)*( BC 2,27 ,K)+BC1,26,K)+B( 1,27 ,K+1) BC! . 2/7 


ze . 5%*(A(6,NB+8)+A(7,NB+9))-5. *B(1,27,K)) 
BOUT( 27, 27,K)=FO(2)*(B( 26,27 ,K)+B(27,26,K)+B(27,27,K+1)+ 
x B(27,27,K-1)+. 5*(A( 16 ,NB+8)+A( 15 ,NB+9) )-5. *B(27,27,K)) 


13 CONTINUE 
* H. THE EXTERIOR CORNERS (COARSE INTERFACE) 


BOUT(1,1,1)=FO(2)*(B(2,1,1)+B(1,2,1)+BOle 1.2) 4 S( ACen. 
* AC7,NB=1) )}4BLEG@2) “DINE =( 456+ EB G2) ob lees) 


BOUT(1,1,10)=FO(2)*(B(2,1,10)+B(1,2,10)+B(1,1,9)+. 5*(A(6,NB)+ 
* AC 7 ,NB-1))+BI(2)*TINF-(4. +BI(2))*B(1,1,10)) 
BOUT( 1,27, 1)=FO( 2)**(B(2,27,1)+B(1,26,1)+B(1,27,2)+. 5*(A(6,NB+8)+ 
x AC 7 ,NB+9) )+BI(2)"TINF-(4. +BI(2))*B(1,27,1)) 
BOUT(1, 27, 10)=FO( 2)*(B(2,27,10)+B(1,26,10)+B(1,27,9)+. 5% 
(A(6,NB+8)+A(7 ,NB+9))+BI(2)*TINF-(4. +BI(2))#B(1,27,10)) 
“BOUT(27, 1,1)=FO(2)*(B(26,1,1)+B(27 52, 1)+B(27,1,2)+. 5%*(A( 16 ,NB)+ 
ss AC 15 ,NB-1))+BI(2)*TINF-(4. +BI(2))*B(27,1,1)) 
Bour(27, 1,10)=FO(2)*(B(26,1,10)+B(27,2,10)+B(27,1,9)+. 5* 
(A( 16 ,NB)+A(15,NB-1))+BI(2)*TINF-(4. +BI(2))*B(27,1,10)) 
‘BOUT( 27, 27, 1)=FO(2)*(B(26,27,1)+B(27 ,26,1)+B( 27 ,27,2)+. 5% 


*e (A( 16, NB+8)+A(15, NB+9))+BI(2)* TINF-(4. +BI(2))*B(27, 27,1)) 
BOUMG@Z7.275 10)=FO( 2): ‘(B( 26,27 ,10)+B(27,26,10)+B(27,27,9)+. 5* 
ye (AC 16, NB+8)+A(15, NB+9) )+BI(2)* ‘TINF-(4. +BI(2))* <B(27, 27, 1008 


* I. THE BOTTOM OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 20 I=10,18 
DO 20 J=NC,NC+8 
TC=0. 
DO 21 IC=3*(1-10)+1,3*(1I-10)+3 
DO 21 JC=3*(J-NC)+1,3*(J-NC)+3 
Dal TC=TC+C(IC,JC,8) 
BOUT(I,J,4)=FO(2)*( BC I+1,J,4)+B(1-1,3,4)+B(1,J+1,4)+B(1,J-154) 
ze +B(1,J,5)+TC/6. -6. 5*B(1I,J,4)) 
20 CONTINUE 


* J. THE SIDE OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 22 J=NC,NC+8 
DO 22 K=2,3 
TC1=0. 
TC2=0. 
DO 23 JC=3*(J-NC)+1,3*(J-NC)+3 
DO 23 KC=3*(K-2)+3,3*(K-2)+5 
TC1=TC1+C(1,JC ,KC) 
23 «*TO2—ne C27 IG es 
BOUT(9 ,J,K)=FO(2)*(B(8,J,K)+B(9,J+1,K)+B(9,J-1,K)+B(9,J,K+1)+ 
zc B(9,J,K-1)+TC1/6. -6. 5*B(9,J,K)) 
BOUT( 19 ,J,K)=FO(2)*( B(20,J,K)+B(19 , J+1,K)+B(19,J-1,K)+ 
* B(19,J,K+1)+B( 19 ,J,K-1)+TC2/6. -6. 5*B(19,J,K)) 
22 CONTINUE 


* K. THE ENDS OF THE EXCLUSION ZONE (FINE INTERFACE) 
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ale 
re 


DO 24 I=10,18 
DO 24 K=2,3 
TC1=0. 
TC2=0. 
DO 25 [0=3*(1-10)+1,3*(1-10)+3 
DO 25 KC=3*(K-2)+3,3*(K-2)+5 
Tel—=M1+C( 1G. 1 Ke) 
2) NG2Z—TG2Z+6C0C.27 ,KC) 
BOUTCI,NC-1,K)=FO(2)*( BCI+1,NC-1,K)+BC1I-1,NC-1,K)+BC1I,NC-2,K) 
* +BCI,NC-1,K+1)+B(1,NC-1,K-1)+TC1/6. -6. 5*B(I,NC-1,K)) 
BOUT(C I ,NC+9 ,K)=FO(2)*( BC I+1,NC+9 ,K)+BC I-1,NC+9 ,K)+BC I ,NC+10 ,K) 
* toc l NCro Kt +B Ll, NG+oeK=1)+9C2/6. <6. 5*BC1 ,NC+9.K))) 
24 CONTINUE 


L. THE SIDE EDGE OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 26 J=NC,NC+8 
TC1=0. 
TC2=0. 
DO 27 JC=3*( J-NC)+1,3*(J-NC)+3 
Met—Teltc( 1, JC.1)+2*6( 1.96, 2) 
Pye TC2=TC2+C( 27 , JC, 1)+2"C(27, JC, 2) 
BOUT(9,J,1)=FO(2)*(B(8,J,1)+B(9,J+1,1)+B(9,J-1,1)+B(9,J,2)+ 


* ely Ort pl Zenit SS p pie 2) yeBCo..) 21) 
BO enon. HOC 2 CBC 200) eB erle  )rB Coe =1 1) BC19 2 )+ 
* eZ O2th i 2) hUNR—( 5S. oth 2) B19.) 51) ) 


26 CONTINUE 
M. THE END EDGE OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 28 I=10,18 
mS1=0). 
TC2=0. 
DO 29 IC=3*(I-10)+1,3*(1-10)+3 
MGil—-16 IHCC1G. 1, 1)+2*C010 , 1,2) 
Pome) TC2=1C24+C(1C,27,1)+2*C(1C,27,2) 
BOUT( I ,NC-1,1)=FO( 2)*(B(I+1,NC-1,1)+B(I-1,NC-1,1)+B(1I,NC-2,1) 
x +B(I,NC-1,2)+TC1/6. +BI(2)*TINF-(5. 5+B1I(2))*B(I,NC-1,1)) 
BOUT( I ,NC+9 , 1)=FO(2)**(B( I+1,NC+9,1)+B( I-1,NC+9,1)+B(I,NC+10, 1) 
+B(I,NC+9 ,2)+TC2/6.+BI(2)*TINF-(5. 5+BI(2))**B(I,NC+9,1)) 
28 CONTINUE 


RETURN 
END 


SOURCE LISTING OF SUBROUTINE COR FORTRAN 
SUBROUTINE COR(A,B,AOUT) 
REAL A,B,AOUT 
DIMENSION A(21,36),AOUT(21, 36),B(27,27,10),FO(3),BI(4) 
COMMON NB,NC,YARC,SUM,TINF,BI,FO 


A. THE CORNERS 


ay, 


AOUT( 1.17)=0: 
AOUT( 1,36)=0. 
AOUT( 21, 1)=0. 
AOUT( 21, 36)=0. 


*" Bo TRESENDS 


DO ioe 
AOUT( I, 1)=FO(1)*(AC I+1, 1)+AC1I-1,1)+AC1,2)+TINF*BI(1)-(3+BI(1))* 
5 (Cl) 
AOUT(1,36)=FO(1)*(A(I+1, 36)+A(I-1, 36)+A(I, 35)+TINF*BI(1)-(3+BI(1) 
eG 30) 
1 CONTINUE 


aw. 


s C. THE SIvES 


DO 2 J=2 55 
AOUT(1, J)=FO(1)* (AC1,J+1)4+AC1,J-1)+AC2,J)+TINF*BIC1)-( 33ers 
vs A( ica 
AOUT( 21, poeta *(AC21,J+1)4AC 21, J-1)4AC 20 (J) 4T INE bie ee 
¥ BiCiIgpAC2i- I) 
2 CONTINUE 


* D. THE MID POINTS 


NMIN=NB-1 
NMAX=NB+9 
DO 3 I=2,20 
DO 3 J=2,35 
IF (J. GE. NMIN. AND. J. LE. NMAX. AND. I. GE. 7. AND. I. LE. 15) GOTO 3 
IF (J. GE. NB. AND. J. LE. (NB+8). AND. (I. EQ. 6. OR. I. EQ. 16)) GOTO 3 
AOUT(I ,J)=FO( 1)*( AC I+1,J)+A(1I-1,J)+A(1,J+1)+A(1,J-1)+BI(1)*TINF 
a =(4+BEGE))-ACi® 
3 CONTINUE 


%* E. THE MEDIUM ENDS TRANSITION 


DO 4 I=7,15 
TA=0. 
TB=0. 
DO 5 IB=(1-7)*3+1,(1-6)*3 
TA=B(IB,1,1)/2.+TA 
TA=B(1B,1,10)/2.+TA 
TB=B(IB,27,1)/2.+TB 
TB=B(IB,27,10)/2.+TB 
DO 5 K=2,9 
TA=B(IB,1,K)4TA 
5  TB=B(IB,27,K)+TB 
AOUT( I ,NMIN)=FO( 1)*( AC I+1,NMIN)+A(I-1,NMIN)+A(1,NMIN-1)+TA/18. + 


zc BI(1)*TINF-(4. 5+BI(1))*A(I,NMIN)) 
AOUT( I ,NMAX)=FO( 1)**(A( I+1,NMAX)+A(I-1,NMAX)+A(I ,NMAX+1)+TB/18. + 

Xe BI(1)*TINF-(4. 5+BI( 1) )*A(I,NMAX) ) 

4 CONTINUE 


* F. THE MEDIUM SIDE TRANSITION 
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DO 6 J=NB,NB+8 

TA=0. 

TB=0. 

DO 7 JB=(J-NB)*3+1,(J-NMIN)*3 
Sp aya gy" 
TA=B(1,JB,10)/2.+TA 
TB=B(27,JB,1)/2.+TB 
TB=B(27,JB,10)/2.+TB 
BO07k=259 

TA=B(1,JB,K)+TA 
7 TB=B(27,JB,K)+TB 
AOUT(6,J)=FO(1)*(A(5,J)+A(6,5-1)+A(6,J+1)+TA/18. + 


* BiGh)*LINE=(43 57801) )2AC6, J )) 
AOUT(16,J)=FO( 1)*(AC17,J)+AC16,J-1)+A(16,J+1)+TB/18. + 
* BIC 1)*TINF-(4. 5+BI(1))*A( Oy ))) 
6 CONTINUE 
RETURN 


END 
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APPENDIX B. STABILITY OF RUNGE-KUTTA 
Applying fourth order Runge-Kutta to the semi-discrete, explicit finite difference 
equations has an effect on their stability. The explicit technique has a three dimensional 
stability criterion of (1-6Fo) > 0, or 


Fos (B.1) 


ae 
6 
For prescribed values of Ax and o, this criterion determines the upper bound on At. The 
stability criterion was determined by collecting all of the 7?,, terms, to obtain the form 
of the coefficient. 

A similar expression can be determined by gathering the terms of 7?,, after the 
completion of Runge-Kutta. This requires a good deal of algebra since a fourth order 


polynomial is the result. The general form of Runge-Kutta 1s shown in equation B.2. 


Al =f(T") 
Kr=f(T' +4) 
; n, K2 
KOE lem) (B.2) 
K4 = f(T" + K3) 
Be an (K1 + 2(K2 + K3) + K4) 


6 


the ayes 


subscript will be dropped and the different terms will be indicated by a relative subscript, 


Because of the amount of notation required to derive the coefficient of 7%,, , 
where 000 is the same as i,j,k and 100 is i+ 1,j,k, and so forth. Substituting the finite 


difference equation for f into equation B.2 and using this notation gives 


Koo = Fo(T 90 7 T_yo9 + To10 + To-10 + Toor + Too-1 — ©To00) 


A200 = Alooo +22 —(Kigg + K1_i99 + Koyo + Aloiig + Aloo: + Aloo_; — 6X1 oo) a 


s 
K3q99 = Aloog + =- 5 (K2j99 + K2_199 + K2q19 + K29_39 + K2001 + K299-1 — OK2o00) 


Koo9 = Klooo + FO(K3 409 + K3_199 + K3019 + K39~19 + K3001 + K3q0-1 — 6K 3000) 


A difference between the form of equation B.2 and equation B.3 is that, in this case, one 
can take advantage of f for a finite difference equation being a linear operator. This 
property of f(7 + KA) = f(7) + f(A) allows greatly simplifving the algebra. 

Now all that is required is to group the terms of 7,,, and find the coefficient. The 


simplest is the KI term, which from equation B.3 is seen to have a coefficient of 


K logo Bas Pac 6Fo 


. (B.4) 
Kl yoo = lFo 


The Aly, 1s the case general. Any of the other five adjacent differences, -100, 010, 0-10, 
OO] and 00-1, will have the same coefficient of 74, as 100. This simplifies finding the 
coefficients for the rest of the terms. 

For Kk2 it becomes a little more difficult. The contribution of K1,,) is already known. 


The values for K1j, tvpe terms 1s given in equation B.4, which yields, 


K2q99 + — 6Fo + 52 (Fo + Fo + Fo + Fo + Fo + Fo — 6( - 6Fo)) ae 


ae Ohooe ho. 


The next term to expand is K3. This will be done in several parts, the 


Ko and A2,., terms are already known. Equation B.6 shows the expansion for a K2hq, 


term. 
ae Fe Da F y Dad Da 
K2 1099 = Al yoo a (Alaoo+ Alooot Aliyot Al yiyot Alioi+ A] j9-,— 641 190) 
| { | { | | | (B.6) 
Fo QO —6fo 0 0 @) 0 — 6Fo 


Kee Eo OF o- 


Since there are six terms of the form K2,, the total coefficient for K3 1s, 


K3o99 7 — 6Fo + #2 (6(Fo — 6Fo") — 6( — 6Fo + 21Fo’)) (B.7) 


SS =6Fo 2 lnor ellno: 


The final term is K4, which has a long expansion. Again it will be broken into parts. 
The Alo and A3,. terms are already known from the above. It 1s necessary to expand 


ime A>, terms, as shown below, 


103 


LO d 
K3109 = Alyoot+ 5- 5 (K2a99+ K2oo9+ A2i 9+ K2y 9+ K219)+ K2)9_1— 6K21q9) 


B.8 
J II III IV IV IV IV Vv ea) 


Equation B.8 has five types of terms, as identified by the Roman numerals. The type J, 
type III and type V terms are already defined by equations B.4 and B.7. Expanding the 


type II term gives, 


ip ; , 
2399 = Algoot 7 (Alg00+ Alioot Alaiot Ala—io+ Algo: + K1g9-1— 641 290) 


! Ue ! | | ! { 
0 0 Fo QO QO 0 0 0 
, 0° 
K2209 > 


Expanding the tvpe IV term gives, 


SCAG a K1y 9+ 2 5 (K1qi0+ Kloot Aljaot+ Alioot Aliyy+ Al,)~,— 641) 10) 


{ { { 1 { { { 1  (B.10) 
0 0 Fo 0 Fo 0 0 0 
on Fos 


Combining equations B.4. B.7, B.8, B.9 and B.10 gives, 


fe, Fe 


KSio — Home — 6Fo + 21 Fo* + 4(Fo*) — 6(Fo — 6Fo’) ] 


B.\) 
123Fo° 


Sho = 6fo- = 7 


The coefficient for the terms in K4 of equation B.3 can now be found by using equations 
B.4, B.7 and B.11. This gives, 


3 
Koo + — 6Fo + Fol 6(Fo — 6Fo? + +#=£2—) — 6( — 6Fo + 21Fo? — 81Fo°)] 
(B.12) 
=, GEE Aone ast ae 


The final step is to gather all of the terms from equation B.2 into the Runge-Kutta 
sum of equation B.]. This results in the coefficient of 77,, , and is given in equation 
Bre: 
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-+ 1 ~ 6Fo + 21Fo? — $4Fo? + 34 Fo" (B.13) 
This fourth order polynomial must satisfy two conditions for stability, the traditional 
explicit, that it be greater than zero and an additional condition, that it be less than one. 
The traditional lint prevents oscillations about a solution which violates the laws of 
thermodynamics. The second limit is based on the law of heat conduction, for if the 
coefficient is greater than one, it is possible for the temperature at a point to rise, even 
though every other point is at a lower temperature. 

The above limits the maximum value of the Fourier number for a three dimensional 
problem to 0.37, which is better than the normal value of one sixth by over a factor of 
two. While not a limit, operating above the minimum point of the Runge-Kutta curve, 
Which occurs at a Fourier number of about 0.22, causes a rapid loss of accuracy. This 
should be apparent, for at this point, taking a larger time step causes smaller temper- 
ature changes instead of the expected larger. A plot of the temperature coefficient versus 
Fo is shown in Figure 36, where it is interesting to note that, for Runge-Kutta, the co- 
efficient never drops below zero. This implies that for a three dimensional problem using 
Runge-Kutta, taking too large a time step will not cause solution oscillation, but rapid 
increase in error. Also shown for comparison is the coefficient for the implicit method. 
Even though this method 1s unconditionally stable, it does not have greater accuracy for 


the Fourier numbers where the other two methods are stable. 


SPEER MEG: - 2 Fo’ +> Fos (B.14) 
eee l0Fo? — 22 Fo’ +42 Fo! (B.15) 


Similar calculations may also be performed for one and two dimensions. These re- 
sults are shown in equations B.14 and B.15, respectively. The comparisions of the tem- 
perature coefficients versus Fourier number are shown in Figure 37 for the one 
dimensional and Figure 38 for the two dimensional finite difference equations. The use 
of Runge-Kutta improves the stability of the explicit finite difference in all cases, how- 


ever the improvement is greater the higher the dimensions. 


THREE DIMENSIONAL 


=) 
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LEGEND 
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Figure 36. Temperature Coefficient Versus Fourier Number for Three Dimesions 
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ONE DIMENSIONAL 


LEGEND 
EXPLICIT 
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Figure 37. Temperature Coefficient Versus Fourier Number for One Dimension 
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TWO DIMENSIONAL 
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Figure 38. Temperature Coefficient Versus Fourier Number for Two Dimensions 
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APPENDIX C. SOURCE LISTING AND EXPLANATION OF AUXILIARY 
COMPUTER PROGRAMS 


A. PLOTTING PROGRAMS 

There is one major plotting program, named PLOT, and several variants which have 
been modified to handle particular cases. The source listing that follows is the most 
general. This program is designed to read the binary output files produced by the weld 
modeling programs and provide either contour or three dimensional drawings of the 
temperatures along a plane slicing the weld area. Since this program grew as the need 
arose, it 1s not well organized. However, it is an interactive program and can be easilv 
used. 

There are three basic presentations available. A temperature contour at a horizontal 
surface at anv depth, a temperature three dimensional plot at a horizontal surface at any 
depth and a contour of a vertical surface at anv X location. A secondary option 1s to 
draw temperaiure profiles for comparison, which was onlv used when validating the weld 
model. 

1. User Guide 

To use this program, vou first must have available the output files from running 
the program WELD, WELDLF or WELDMA. Insure that vou have at least 1500K of 
storage by typing GETSTOR 1500K in CMS. Your profile exec should include the 
NONIMSL library, since one routine is called from this library. Then type DISSPLA 
PLOT, to execute this DISSPLA program. After tvping enter in response to the ques- 
uon about file definitions, the program will be loaded. 

There are two ways to displav the output of the program, the first 1s onto the 
screen of a tektronics and the second Is to create a metafile and then displav the output 
using DISSPOP. The later is done by choosing the COMPRS opuon, which then allows 
using the laser printer for higher quality graphics. The next question has vou select 3D 
plots or profiles, in general vou will select 3D plots. You will next be asked the contour 
interval, 500 is a good tvpical value. At this point you are in the program loop so you 
can try another value later. The last preparatory question allows vou to add a specific 
label to the graph. This label can be 60 characters long and should be ended with a ‘S’ 


or else the title will not be centered. 
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The next choice is which data you desire to see. The first three choices are taken 
from files which saved a specific surface at half second intervals. The last three choices 
are from the final data saved and hold the entire temperature matrices and either hori- 
zontal X-Y planes or vertical X-Z planes may be looked at. If you select a running 
surface, you will be asked for a time. Since there is one surface saved every 0.5 seconds 
of simulation, you can ask to see any of these. If you attempt to look at a surface that 
hasn't been saved vet, then the program will fail. For the last three choices, you will 
have to specify the type of plot desired and the location of the plane to plot. As indi- 
cated in the question, these locations are always referenced to the nodes of the surface 
vou are looking at. If vou are looking at an X-Y surface you will be asked if vou want 
a 3D plot. A 3D plot draws the temperature as a third dimension over the surface, 1n- 
stead of a contour map. 

If vou have selected Tektronics, the plot will now be displayed. You can then 
choose to draw another plot. If vou have selected compress, a metafile is created. After 
vou have made all of the plots vou desire, vou can use DISSPOP to send vour plots any 
Where vou choose. Just type DISSPOP and enter and then answer the prompts. 

2. Program Description 

The program uses the package of graphics programs called DISSPLA [Ref. 20] 
and this discussion assumes that the reader 1s familiar with these routines. The first 
section of the program defines the titles that are placed on the head of the graphs. By 
use of an equivalence statement and internal formatted writes these titles are modified 
as needed to support correctly identifying each graph. 

The second portion of the program inputs the major options selected by the 
user. The first choice is to use a tektronics or to use DISSPOP, the later being run after 
this program, using the metafiles created by the COMPRS subroutine. The second 
major Opuon is to create three dimensional temperature profiles or to compare line 
profiles. This later choice is a simple routine which makes use of the NONIMSL routine 
PLOTD. For normal plotting the next choice will be the contour interval. The contour 
interval should be large enough to keep three dimensional graphs on scale and contour 
plots to not have an excessive number of lines. Experience has shown that 500 degrees 
is a good value. The next prompt allows inputting a fourth line of title to include addi- 
tional information about this picture, such as arc power, speed, etc. The last question 
of the input section asks which type of plot is desired. 

Since the Fine and Medium zone horizontal surfaces are both square, the pro- 


gram 1s Written with them using the same calls. Since the scale is different, the first step 
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is to set the scale variable. This is only used if a Fine or Medium horizontal plot is being 
made. 

The next section is a large IF-THEN-ELSE statement which handles the two 
major cases, running output or final output, for graphing. The first half is running 
Output, so the user 1s asked at what time the output is desired, a time which will alwavs 
be truncated to the nearest 0.5 interval of time. If either a Fine or Medium surface 1s 
to be read, then the file SURF 1s read until the appropriate surface is inputted. If the 
X-Z profile is to be read, then the file CUT is read until the appropriate vertical surface 
is inputted. In this last case, the matrix is read in reverse order for the Z direction to 
allow proper orientation on output by the contour routine. The variable [MAT is used 
to identify the type of output geometry to use later in the program. At this point the 
AREA2D DISSPLA subroutine is called to define the plot area. 

The second half of the IF-THEN-ELSE statement handles the three cases where 
the FINAL output will be used to create the plots. The three dimensional arrays (two 
dimensional for the Coarse zone) are read in. If the Coarse zone is chosen. there is onlv 
One surface since it is two dimensional, and the AREA2D 1s called and IMAT set. For 
the Fine and Medium Zones, either the horizontal X-Y plane or vertical X-Z plane must 
be selected for plotting. Next, the location of the plane must be specified, that 1s Z for 
the X-Y plane and Y for the X-Z plane. This value is in nodal position, not physical 
distance. The program will convert nodal position to distance when drawing the plot. 
Again IMAT is set for the type of plot and AREA2D is called to set the plotting area. 

At this point the IF-THEN-ELSE block has ended and the option of plotting 
a three dimensional profile versus a contour plot is given for horizontal Fine or Medium 
surfaces. If a contour plot is chosen the sequence of DISSPLA subroutines 1s called to 
produce the desired plot. The four types of surfaces that can be plotted are determined 
by the variable IMAT, which 1s specified in the variable definitions at the start of the 
program. When the contour plot is done, the user is allowed the option of creating an- 
other plot. If for the Fine or Medium horizontal plots the three dimensional plot was 
chosen, then the DISSPLA routines for this plot are called. Again the program allows 
the option of producing another graph. 

The last section of the program 1s for the special case of producing comparative 
temperature profiles in the fine zone from several different sources. The input files 
should be the Final files from running the Weld programs or the Rosenthal verification 
program. These are binary files and should have different names. Either the X or the 


Y axis mav be chosen for the graph. This same axis will be used for all graphs. It is 
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possible to plot several temperature profiles from the same file, all at different Y-Z (X-Z) 
locations. Upon completion, the program loops to allow producing additional plots. 
a. Fartants 

There are two major variants to the plotting program. There is a modifi- 
cation, PLOTL, for plotting the larger arrays associated with the WELDC program 
which is used for cooling rate calculations. The only difference in this program is that 
the Fine and Medium arravs have been enlarged and the indices changed were appro- 
priate. There is a modification, PLOTYZ, for plotting Y-Z vertical planes vice the 
standard X-Z vertical planes. The onlv difference in this program 1s that for Final Fine 
and Medium plots, the Y-Z plane can be selected vice the X-Z plane. 

3. Source Listing of Program PLOT 


PROGRAM PLOT 


Yetervededevevedesevedevededevedovedededededededesevedevescslievedevescscdeveukevedesecdededesedcadevedesdcveskcdevedeveseutvesesedevedescstcscvevesieae 


* THIS PROGRAM PLOTS CONTOURS OF THE OUTPUT DATA FROM RUNNING THE * 


* PROGRAM WELD. THE FOLLOWING VARIABLES ARE DEFINED: * 
seats >: THE TEMPERATURE DIFFERENCE BETWEEN CONTOURS * 
* Ni SPECIFIES Witch PIO iE AvTNG 5 Pinko? * 
a2 : SPECIFIES WHICH PLOT HEADING IS SECOND * 
* ON > LIST OF LENGTHS OF THE HEADINGS * 
* HEAD : THE LIST OF PLOT HEADINGS * 
* HD,DEPTH, TIMES ,TIME8 ,H9Y,HIOY ARE ALL USED TO MODIFY HEADINGS * 
*%* IMAT : INTEGER SPECIFYING PLOT TYPE BEING OUTPUT (DIMENSIONS) * 
* VALUE PLOde® TY PE ASSOCIATED MATRIX % 
* i KY PLOT FOR SRINE OR MEDIUE ZMAT * 
* Z 5G Ol VO rink XZF ¥ 
* 5 XY PLOT FOR COARSE ZCMAT 3 
* 4 42 PLOT FOR MEV XZM * 


dedvdevedededcsevdedeucvesevedeseveveakedeskevevescotscalkclcscvevevencotcveskeveveskacseacaescseseakcakeakeakaleskakdevcsesle desks sede sleveaevleslcakeae seve 


DIMENSION ZMAT(27,27),XZF(27,8), ZCMAT( 21,36) ,XZM(27,10),N(10) 
*.C(27,27,8) ,BC27 , 210) , MeFi? Aes) ONC AAG, 

CHARACTER HEAD(12)*50,TIME4**4, TIME8*4 , DEPTH*"2 , H9Y**3 ,H10Y*3,H(50) 
** HD*50,ANS , TEMPS , USERHD**60 , AXS , FNAME*'8 , XLABEL**60 , YLABEL*60, 
*LGDTXT*20 

EQUIVALENCE (HD,H(1)),(H(8) ,TIME4) ,(H(23),TIME8) ,(H(9) , DEPTH) 

* (H(40) ,H9Y) ,(H(42) ,H10Y) ,(H(37) , TEMP) 
COMMON WRK( 20000) 
DATA N/29,31,33, 19,21, 23,23, 34 


45,47/,HEAD/ 
*'FINE ZONE SURFACE TEMPERATURE 


' "MEDIUM ZONE SURFACE TEMPERATURE' , 


*' PROFILE OF X-Z TEMPERATURE AT ARC','TIME =  . SECONDS', 
*'FINE ZONE TEMPERATURE' ,'MEDIUM ZONE TEMPERATURE’ , 

*' COARSE ZONE TEMPERATURE’ ,'DEPTH = MM TIME =  .  SECONDS', 
*'FINE PROFILE OF X-Z TEMPERATURE AT Y = MM', 

*'MEDIUM PROFILE OF X-Z TEMPERATURE AT Y = MM', 
*'TEMPERATURE BETWEEN CONTOUR LINES = DEG han 


*'FINE ZONE TEMPERATURE PROFILESS' / 


C OFFER OPTIONS OF DISPLAY TO USER 


Oe 


WRITEC (1X. A)" ) 
“i THIS PROGRAM DISPLAYS THE OUTPUT FROM THE PROGRAM WELD. ', 
*'THE OUTPUT CAN BE DIRECTED TO A TEK618 PLOTTER OR TO A COMPRS', 
*'FILE FOR LATER OUTPUT. THE OUTPUT ARE CONTOUR PLOTS OF THE', 
*'FINE SURFACE, THE MEDIUM SURFACE OR THE XZ PROFILE AT THE ARC’, 
*'AT .5 SECOND INTERVALS OR ANY POSSIBLE HORIZONTAL OR XZ SLICE’, 
*'OF THE FINAL TEMPERATURE DISTRIBUTION. FOR FINE AND MEDIUM’, 
*' HORIZONTAL SURFACES IT IS POSSIBLE TO DRAW 3D PLOTS. FOR THE', 
*'FINE ZONE ONLY IT IS POSSIBLE TO DRAW COMPARATIVE TEMPERATURE' , 
PROFILES.’ ,' ', 
*' PLEASE SELECT PLOTTING DEVICE: (1) TEXTRONIX 618 (2) COMPRS' 

READ (*,%*) IPLOT 

IF (IPLOT. EQ.1) CALL TEK618 

IF (IPLOT. NE. 1) CALL COMPRS 
1 WRITE(*,*)' DO YOU WANT: (1) 3D OR (2) PROFILES?' 

READ(*,*) L 

IF (L.EQ.2) GOTO 1000 

WRITE(*,'(/,1X,A)') 'WHAT IS THE CONTOUR INTERVAL? ' 

READ(*,*) T 

HD=HEAD( 11) 

WRITE(TEMP,'(F5.1)') T 

iw / ce 

HEAD( 11)=HD 

WRITE(*,*)' TYPE IN ANY ADDITIONAL GRAPH LABEL DESIRED, TERMINATE' 

WRITE(*,*)' WITH A $. (60 CHARACTERS MAX) ' 

READ(*, '(A60)' )USERHD 

WRitE(@= (1X ,A) ) °., 

*'SELECT ONE OF THE FOLLOWING OPTIONS: ', 

*'(1) FINE SURFACE AT .5 SECOND INTERVAL’ , 

*'(2) MEDIUM SURFACE AT .5 SECOND INTERVAL’, 

*'(3) X-Z PROFILE OF ARC AT .5 SECOND INTERVAL’ , 

*'(4) FINAL FINE ZONE TEMPERATURE’ , 

*'(5) FINAL MEDIUM ZONE TEMPERATURE’, 

*'(6) FINAL COARSE ZONE TEMPERATURE ' 

READ(* ,**) IPLOT 

IF (IPLOT. EQ. 1. OR. IPLOT. EQ. 4) THEN 

SCALE=1- 
ELSE 
SCALE=3. 
ENDIF 


C INPUT THE DATA FOR THE OPTION SELECTED 


IF (IPLOT. LE. 3) THEN 
N1=IPLOT 
N2=4 
WRITE(*,*) ' WHAT IS THE DESIRED TIME?' 
READ(*,*) TIME 
TIME=( INT( 2*TIME) /2. ) 
HD=HEAD(4) 
WRITE(TIME4, '(F4.1)')TIME 
ISKIP=TIME/. 5-1 
IF (IPLOT. LE. 2) THEN 
C SURFACE PLOT 
OPEN(1,FILE='SURF' ,STATUS=' OLD' , FORM='UNFORMATTED' ) 
C SKIP OVER EARLIER TIME DATA IN FILE 
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DO 10 I=1,ISKIP 
READ( 1) ZMAT 
10 READ( 1) ZMAT 
READ(1) ZMAT 
IF (IPLOT. EQ. 2) READ(1)ZMAT 
CALL AREAZD(U 77. ) 
IMAT=1 
ELSE 
C FINE XZ PROFILE AT SPECIFIED TIME 
OPEN(1,FILE='CUT' ,STATUS='OLD' , FORM='UNFORMATTED' ) 
BO tel— ie TSKIP+ 
11  READ(1) XZFH 
DO 121 11 
DO 111 K=1,8 
ial XZFC 1 k= Zee ok) 
IMAT=2 
CALL AREA2D@7 alge 2as) 
ENDIF 
ELSE 
C PLOT CONTOURS FROM DATA AT FINISH (FINAL) 
OPEN(1,FILE='FINAL’' ,STATUS='OLD' , FORM='UNFORMATTED' ) 
READ(1) TIME 
READ Gime 
READ(1) B 
IF (IPLOT. EQ.6) THEN 
C PLOT COARSE DATA, SKIP OVER FINE AND MEDIUM DATA IN FILE 
READ( 1) ZCMAT 
HD=HEAD( 8) 
WRITE( DEPTH, '(12)')9 
WRITE( TIME8, '(F4.1)')TIME 
N1=7 
N2=8 
IMAT=3 
CATE VAREA2DCGadew 2) 
ELSE 
C SELECT PROFILE TYPE FOR FINAL DATA PLOTTING 
WRITE(*,*)' SELECT OPTION: (1) HORIZONTAL (2) X-Z PROFILE’ 
READ(* ,**) IOPT 
IF (IOPT. EQ. 2) THEN 
Le WRITE(*,*)' AT WHAT Y LOCATION (1 TO 27)?' 
READ(*,*) J 
IECJ. LT: 1. OR. J.GT.27 COhGms 
ENDIF 
IF( IPLOT. EQ. 4) THEN 
C FINE FINAL DATA PLOT, HORIZONTAL SURFACE 
IFC IOPT. EQ. 1) THEN 
WRITE(*,*)' AT WHAT DEPTH (1 TO 8)?’ 
READ(* ,%) K 
DOA ie 
Doe —e 
14 ZMATC Te CC ae 
IMAT=1 
N1=5 
N2=8 
HD=HEAD( 8) 
WRITE( DEPTH, '(12)')K-1 
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WRITE(TIME8, '(F4. 1)' )TIME 
CALE AREA 77 ) 
ELSE 
C FINE FINAL PLOT, XZ PROFILE 
ne) Ss il ey, 
Oris k—is 
15 XZF(1,9-K)=C(1I,J,K) 
IMAT=2 
N1=9 
N2=4 
HD=HEAD( 9) 
WRITE(H9Y,'(13)')J 
HEAD( 9 )=HD 
HD=HEAD(4) 
WRITE( TIME4,'(F4.1)')TIME 
CAME ARBAZDC 7: . 2.4) 
ENDIF 
ELSE 
C MEDIUM FINAL PLOT, HORIZONTAL SURFACE 
IF(IOPT. EQ. 1) THEN 
WRITE(*,*)' AT WHAT DEPTH (1 TO 10)?' 
READ(*,*) K 
DO 16 I=1,27 
DOmlGs J=1 27 
16 Zo BC ed) 
maT = 1 
N1=6 
N2=8 
HD=HEAD( 8) 
WRITE( DEPTH, '(I2)')K*3-3 
WRITE(TIME8,'(F4.1)')TIME 
CALL AREA2D(7,7) 
ELSE 
C MEDIUM FINAL PLOT, XZ PROFILE 
Dom 7 k= id 
DO) eS ey, 
Wy P77 ieiek )—pC 1. J,11-K) 
IMAT=4 
N1=10 
N2=4 
HD=HEAD( 10) 
WRITE(H10Y, '(13)' )J*3-2 
HEAD( 10)=HD 
HD=HEAD(4) 
WRITE(TIME4,'(F4.1)')TIME 
CALL AREA2D(7. 1,3) 
ENDIF 
ENDIF 
ENDIF 
ENDIF 


G ALLOW OPTION OF A 3D PLOT 
IF (IMAT. EQ. 1) THEN 


WRITE(*,*)' DO YOU WISH A 3D PLOT?’ 
READ(*,'(A1)') ANS 


1 


oa 


IF (ANS. EQ. 'Y') GOTO 100 
ENDIF 


TEMPORARY FIX FOR CONTOUR PLOTTING 
[=10 
IF I=] 
CALL THKFRM(. 01) 


SET UP AND PLOT CONTOUR 


CALL BCOMON( 20000) 
CALL XNAME('X-AXIS DISTANCE IN MM' ,22) 
CALL HEADIN(HEAD(N1),N(N1),2,4) 
CALL HEADIN(HD,N(N2),2,4) 
CALL HEADIN(HEAD(11),47,2,4) 
CALL HEADIN(USERHD,100,2,4) 
IF (IMAT. EQ. 1) THEN 
CALL YNAME('Y-AXIS DISTANCE IN MM' ,22) 
CALL GRAF(0, 3*SCALE ,27*SCALE,0, 3*SCALE, 27*SCALE) 
CALL JGR ID GLO) 
CALI CONMARCZMAT. 27 927-15) 
ENDIF 
IF (IMAT. EQ. 2) THEN 
CALL YNAME('Z-AXIS DISTANCE IN MM',22) 
CALUMGR AmGOnsn 27 <7 620) 
CALL GRID(0,0) 
CAMMGONMARCXZF. 27°87) 
ENDIF 
IF (IMAT. EQ. 3) THEN 
CALL YNAME('Y-AXIS DISTANCE IN MM’ ,22) 
CAI CRAP(O+2)7 -180"0) 27 lon 
CALL GRID(0,0) 
CALL CONMAK( ZCMAT.21.°36.T) 
ENDIF 
IF (IMAT. EQ.4) THEN 
CALL YNAME('Z-AXIS DISTANCE IN MM',22) 
CALLM@GRAF( 0 Gest es 0) 
CALEVGRIDCO,0) 
CATIA CONMAR( KZMee7 = 1057) 
ENDIF 
CALL CONANG( 30) 
CALL CONMIN(4) 
CALL CONDIG(0) 
CALL CONTHN(. 05) 
CALL CONLIN(O, SOLID’ ; “HARES 
CALERCONEEN( 1, DOIk NOLABERS. laa) 
CALL CONLIN(2,'DASH' , 'NOLABELS' ,1,3) 
CALL CONLIN(3, ‘DOT’ , 'NOLABELS’ ,1,1) 
CALIY RASRENG= 25) 
CALL CONTUR(4, ‘LABELS’ , ' DRAW’ ) 
CALL ENDPL(0) 
CLOSE(1) 
WRITE(*,%*) ' DO YOU WISH TO TRY ANOTHER CURVE?’ 
READ(*,'(Al)' )ANS 
IF (CANS=EO) ¥° ) GOTO 1 
CALL DONEPL 
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SOP 


C 3D PLOT INSTEAD OF CONTOUR PLOT 
100 CALL HEADIN( HEAD(N1),N(N1),2,3) 
CALL HEADIN(HD,N(N2),2,3) 
CALL HEADIN( USERHD, 100, 2,3) 
CALL X3NAME('X-AXIS DISTANCE IN MM',22) 
CALL Y3NAME('Y-AXIS DISTANCE IN MM‘ ,22) 
CALL Z3NAME('TEMPERATURE IN DEGREES KELVINS' , 100) 
GALL VOEMSDC 10: , 10. ,8. ) 
CALL VUABS(-15. ,-15. ,12. ) 
CALL GRAF3D(0,3**SCALE, 27**SCALE ,0, 3* SCALE , 27*SCALE, 300,T*4, 
* T*24+300) 
SAGs SURMAT( ZMAT,1,27,1,27 0) 
GOTO 99 
C FINAL FINE ZONE PROFILE PLOTTING 
Mego WRITE(*, (/,1X,A) )' ' 
WRITE(*, (1xX,A)') 


am THIS PROCEDURE ALLOWS PLACING MULTIPLE TEMPERATURE PROFILES’, 
*'ONTO ONE PLOT. THE INPUT FILES SHOULD BE NAMED AS FOLLOWS: ', 
ve! FILE FILENAME A', 


7’ WHERE FILENAME IS WHAT IS NORMALLY CALLED FILETYPE IN VMS. 
ONCE' , 
*'THE AXIS IS SPECIFIED FOR A PLOT IT MAY NOT BE CHANGED, THOUGH', 
*'THE OTHER TWO VARIABLES ARE FREE. YOU WILL BE ASKED FOR A LABEL' 
* "FOR EACH GRAPH.',' ' 
DO 1006 I=1,27 
HOGoexDATA(I)=1-.5 
YLABEL=' TEMPERATURE IN DEGREES KELVINS' 
1001 WRITE(*,*)' WHAT IS THE REFERENCE AXIS, X OR Y?' 
READ(*,'(A1)') AXS 
Ret AxSeNE. Y .AND.AXS.NE. X°) GOTO 1001 
1002 WRITE(*,**)' WHAT IS THE FILENAME FOR THIS CURVE?' 
READ(*, '(A8)')FNAME 
OPEN( 1, F ILE=FNAME , STATUS='OLD' , FORM='UNFORMATTED' ) 
READ( 1) TIME 
READ(1) C 
CLOSE(1) 
WRITE(*,*)' WHAT IS THE LABEL FOR THIS CURVE? (20 CHARACTERS MAX' 
WRITE(*,*)' AND MUST END IN A $)' 
READ(*,'(A20)') LGDTXT 
IF (AXS. EQ. 'Y') THEN 
XLABEL='Y AXIS DISTANCE IN MMS' 
WRITE(* ,**) 
WRITE(*,**)' WHAT ARE THE DESIRED X,Z COORDINATES? (1-27,1-8)?' 
READ (*,*) I,K 
DO 100 J=1, 27 
fee YDATA(J)=C(1,J,K) 


ELSE 
XLABEL='X AXIS DISTANCE IN MMS' 
WRITE(* ,*) 


WRITE(*,*)' WHAT ARE THE DESIRED Y,Z COORDINATES? (1-27,1-8)?' 
READ (*,*) J,K 
DO 1004 1=1,27 
1004 YDATA(I)=C(1I,J,K) 
ENDIF 


ey 


WRITE(*,*)' DO YOU WISH ANOTHER PROFILE ON THIS PLOT?’ 

READ(*, ‘(Al)’ )ANS 

IF( ANS. EQ. 'Y') THEN 

CALL PLOTD( XDATA, YDATA,27,. FALSE. ,' LINLIN' , LGDTXT,HEAD( 12) ,XLABEL, 

YLABEL) 

“cone 1002 

ENDIF 

CALL PLOTD( XDATA, YDATA,27,. TRUE. ,' LINLIN' , LGDTXT,HEAD(12) ,XLABEL, 
* YLABEL) 

GOTO 999 

END 


B. ROSENTHAL VERIFICATION PROGRAM 

This simple program uses the Rosenthal solution for a point source to predict the 
temperatures in lieu of the WELD program. It produces a FINAL file of similar format 
to the FINAL file of the weld program which can then be used to produce Final Fine 
or Medium plots of either the three dimensional or contour plots. In addition these re- 
sults may be compared to the WELD program results by using the profile section of the 
PLOT program. 

1. Source Listing of Program ROSEN 


keene: ROSEN 
Vevlevevesesleslestovloclesieste eve eslesleste oule sles Youle we eve ewe sleste we ale se estete este se eres ewes fentle e's Je Je es etedestewowes s'e Je le slevleclonios! esieslc evesie ses erles foto Je s'e sles'e we oe lowe wos lease s'e 
* THIS PROGRAN SOLVES THE ROSENTHAL SOLUTION FOR A LOW POWER HEAT * 
* SOURCE. THIS IS TO VALIDATE A 3 DIMENSIONAL FINITE DIFFERENCE % 
** MODEL OF WELDING ON A THICK PLATE. THE OUTPUT FILE MIMICS A FINAL * 
% FILE GENERATED BY THE FINITE DIFFERENCE MODEL POR THE FINE VANE cS 
* MEDIUM GRIDS AT TIME = 10 SECONDS. THE FOLLOWING PARAMETERS WERE ~* 
7 USE IWe * 


** ARC VELOCITY = 4& MM/S VOLTAGE = 30 VOLTS CURRENT = 265 AMPS * 
ve EFFICTENGY — a0G5 QARC = 39.75 WATTS i 
* K = 53 W/M-K C = 4.23 E6 J/M3-K FINE GRID SPACING = .001 M * 
* TO = 300 K MEDIUM GRID SPACING = .003 M % 
% ne pee Ss ULE POORER ICS a 


sloalenton! fevlestestos wosleslc steverle sfevlovevloslestovlesios wenloclenlenloclsclactonioctenianlionlonlecloaionleniocfonie a! evedevevotevevolevovevodesoveltetetolet: sleveveveveveveveve 


DIMENSION MiG, 27.8). Bae? 727) 
OPEN(1,FILE='ROSEN' ,STATUS='NEW' , FORM='UNFORMATTED' ) 
WRITE(1) 10. 
DO 1 I=1,27 
DO ee eee 
DO 1 K=1,8 
X=(I-14)**. 001 
Y=(J-17)*. 001 
Z=(K-1)*. 001 
R=SORT( X°X+Y"Y+Z"Z) 
IF (R.EQ. 0.) THEN 
AC(I,J,K)=1000. 
ELSE 
A(I,J,K)=300. +. 1144*EXP( -160. *(Y+R)/R 
ENDIF 
1 CONTINUE 
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WRITE(1)A 
DO 2 I=1,27 
DO 2 J=1,27 
Don 2k=1. i110 
X=(1-14)*. 003 
Y=(J-15)**. 003 
Z=(K-1)*. 003 
R=SQRT( X*X+Y*Y+Z*Z) 
IF (R.EQ.0.) THEN 
B(I,J,K)=1000. 
ELSE 
B(1,J,K)=300. +. 1144*EXP( -160. *(Y+R)/R 
ENDIF 
2 CONTINUE 
WRITE(1)B 
END 


C. MISALIGNMENT PROGRAM 

To study the result of tracking off of a weld seam, the main program WELD and the 
FIN subroutine were modified to simulate a seam. These modified versions are called 
WELDMA and FINMA, with the MA indicating misalignment. The WELDMA pro- 
gram differs in that it can position the arc in the x as well as the y direction. This mo- 
dification is only in the portion of the code which adds the energy from the arc, and is 
included below for reference. The FINMA subroutine has been modified to include a 
new variable, MELT, which determines if the nodes adjacent to the seam has melted or 
more melting has occured, than MELT for that node is set equal to TRUE and the 
heat conduction across the seam is calculated normally. If the seam has not melted. 
than the heat conduction for the nodes adjacent to the seam are recalculated, using zero 
thermal conductivity. This modification is only made to the surface and the interior 
nodes, which is listed below. 

1. Source Listing of WELDMA Modifications 


PROGRAM WELDMA 
*kk Modification of how arc energy is added *** 
C POSITION HEAT SOURCE AND CALCULATE VOLUME WEIGHTING FACTOR SUM!’ 


YARC=YVEL*T IME+7 3 -NB*°9 -NC*3 
IF (TIME. LE. 15.) THEN 
XARC=13. 5 
ELSE 
XARC=13. 5+XVEL*(TIME-15. ) 
ENDIF 
SUM=0. 
DO 1 J=10,23 
DO 1 I=10,27 
SUM=SUM+EXP( -. 10625%*( ( J-YARC) **2+( I-XARC)**2) )%. 5 
DO 1 K=1,3 
SUM=SUN+EXP( =. 106257-( (J-YARC)***2+( 1-XARC)**24+K**2 )) 


M9 


1 CONTINUE 
SUM=SUM/Q 


C ADD THE HEAT EFROMS IE ARG 


DO 2 J=10,23 
DO 2 I=10,27 
C(I,J,1)=CC1,J,1)+EXP( -. 10625*( ( J-YARG)#*2+( 1 -XARC )%42 7 ee 
DO 2 K=2,4 
C(CI,J,K)=C(I,J,K)+EXP(-. 10625*( ( J-YARC )**2+( I-XARC )**2+(K-1) 
ve 2) )/SUM 
2 continue 


2. Source Listing of FINMA Modifications 
wee Modification to track status of the seam *** 
C THE INTERNAL NODES 


DO 1 I=2,26 
DO 1 J=2,26 
DO 1 K=2,7 
COUT(1,J,K)=FK(C(1,J,K) )*B0(3)*(GK(C(141,J,K), CCl, J skp 
*% GK(C(I-1,J,K),C(1,J,K))+GK(C(1,J+1,K) ,CCI,J,K))+GK(C( 1 eee 
* ~§6 OC 1.5K) )+GK(C(1, J, KAT) SOC 1, J,K) )tGKCG(1 J K- 1) Gt ee 
1 CONTINUE 
C RECALCULATE THE DIFFERENCE EQUATION FOR THE SEAM IF NOT MELTED 
DO 10 J=2,26 


DO 10 K=2,7 
IF ((C(13,3,K). GEe 1973. WOR) (COL SI KO mGEateaTC an 
ze MELT(J,K) = . TRUE. 


IF (NEMEC Kk) Coromia 
COUT(13,J,K)=FK(C(13,J,K))*FO(3)*( 
* GK(C(12,J,K),C(13,J,K))+GK(C( 13, J+1,K) ,C(13,J,K))+GK(C( 13 me 
* (C0(13,J,K))+GK(C(13,J,K+1),G(13;J,K))+GKCG( 13g, K- 1) , Camara 
COUT(14,J,K)=FK(C(14,J,K))*FO(3)*(GK(C(15,J,K),C(14,J,K))+ 
* GK(C(14,J#1,K),C(14,J,K) )+GK(C(14,J-1,K), 
* C(14,J,K))+GK(C(14,J,K+1) ,C(14,J,K))+GK(C(14,J,K-1),C(14,J,K) 
10 CONTINUE 


K) 
)) 


C THE TOP SURFACE 


TINF4=TINF**4 


DO 2 I=2,26 
DO 2 J=2,26 
COUT(1,J,1)=FK(C(1,J,1))*#F0(3)*(GKOCCI+ lI. CCL ee 
*  GKCG(E+1,J,1),C(1.J,1))+GK(CC1, +1 sien 1 ene 
*  GK(C(I,J-1,1),CC(I,J,1))+2. *GK(C(I,J,2) ,C(1,J,1)))+ 
*  BIC3)*(TINF4-C(1,J,1)#*4)+B1(4)*(TINF-C(I,J,1)) 
2 CONTINUE 
C RECALCULATE THE DIFFERENCE EQUATION FOR THE SEAM IF NOT MELTED 
DO) 260 J=2026 
IF ((C(13,J,1).GE.1773. ). OR.CCC1I4, J) GEa cee 
“ MELT(J,1) = . TRUE. 


IF (MELT(J,1)) GOTO 20 
COUTC 13 ,J,1)=FK(C(13 J, 1) ) HOC 3 0 
*  GK(C(12,J,1),C0(13,J,1))+GK(C(13,J+1,1),C0(13,J,1))+ 
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* CRCC( 13, 3-11) 96013. Ieee eNCGC 13 J 2) CC 13a. 10S 
* Bi Soe (TING C13 Jetje (4 jC INF-CC13,J,.1)) 
COUTC 14 351 )=FK(C( 14.91) )*FOCs ye (CK CC 15.,,J, 196014 Je 
paeeGK(G( 14, J+1,1),G(i4,3,1))+ 
* GK(C(14,J-1,1),C(14,J3,1))+2. *GK(C(14,J,2),C(14,J,1)))+ 
% BIC3)*(TINF4-C(14,J,1)**4)+BI(4)*(TINF-C(14,J,1)) 

20 CONTINUE 
D. LACK OF FUSION PROGRAMS 

To study flaws (or regions of lack of fusion) the main program WELD and the FIN 
subroutine were modified to allow creating regions of zero thermal diffusion. In addi- 
tion, an additional file, COMP, was created to produce the history of the temperatures 
behind the arc at 0.25 second intervals. This was one line of temperatures on the sur- 
face, at a y location of three millimeters directly behind the arc. The distance of three 
millimeters was insured since the arc was moving at one millimeter every .25 seconds and 
the grid spacing was one millimeter. When speeds other than four millimeters per second 
are used, this time should be modified accordingly. Thus after five seconds a total of 20 
data sets had been taken. Only the modifications made to the program WELD and the 
subroutine FIN are listed below, since thev were not extensive. 

Two auxiliary programs were used to process the data from the file COMP. The 
first, called TABLE, merely converted the binary data in the file COMP to a readable 
data table. The second program, called FLAW, was used to fit the results of analyzing 
the data from the program TABLE. This was written in basic, since it was helpful for 
the program to be easilv altered while trving to fit the data. The resulting fit from this 
program was reported 1n chapter 5. 

When running the program, it was desired to perform multiple runs over one night. 
This was done by use of an exec. A sample exec 1s included below. This exec would 
copy the startup files from the backup disk B and then setup the input stack. The pro- 
gram then executes, and the exec prints out the desired results. The versions of plot 
listed are identical to the PLOT program but each does a predetermined set of plots. 
The modified startup files are discarded and the originals recopied to the A disk in pre- 
paration for running the next problem. 

1. Source Listing of Program WELDLF Modifications 


PROGRAM WELDLF 
THIS VERSION IS USED FOR INSERTING LACK OF FUSION ZONES IN THE 
FINE ZONE. THE OUTPUT FILES ARE MODIFIED IN NAME AND IT USES A 
MODIFIED FINE ZONE SUBROUTINE FINLF. THE SUBROUTINE FINLF IS WHERE 
THE ACTUAL SIZE AND LOCATION OF THE LACK OF FUSION IS CONTROLLED. 
THE POSITION OF THE LOF ZONE IS CONTROLLED WITH THE VARIABLE LOF. 


LG 1s Ge a Sp Lal Sp ial 
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* Added before first executable statement 


COMMON: 7 UNEUT) dil? Shi K2 so 


* Added before entering the main program do loop 


C INPUT LACK OF FUSION ZONE DIMENSIONAL DATA 
READ(C*<))) speek KO JL 
C INITIAL CONDITION VERIFICATION REPORT 
OPEN (11,FILE='VERIF' , STATUS='NEW' ) 
WRITE (11,1234) TIME,OUT,Q,NB,NC,VEL,STEP,CSTEP,BSTEP,N,DELT, 
eA OK eae 
CLOSE(11) 
OPEN(11,FILE='COMP' ,STATUS='NEW' , FORM='UNFORMATTED' ) 
1264 FPORMAT( TIME = °.F6:2./ GUD — (F620. eer, 
+ ' NB = ',13,/' NC = ',13,/ VEL = ',F6.2,/ STEP = Sheu 
+ ' CSTEP = ',13,/' BSTEP = ',13,/' N= ME DELT = "Si6eee 
+ / tee" 13,’ TO (ToweeeZ = 813, see.) eee 
+ ' LONG. ') 


big 


ule 


oS 
¢ 


ve 


Ca? Cl Gy Cy Car) 


Typical modified Runge-Kutta step, calls subroutine FINLF versus FIN. 
All four of the calls to FIN are changed to FINLF 


FIRST RUNGE-KUTTA ITERATION 


CALL COR(A,B,ASUM) 
CALL MED(A,B,CIN,BSUM) 
CALL FINLF(B,CIN,CSUM, LOF) 


Inserted immediately following completion of the Runge-Kutta sum. 


OUTPUT COMPARISON OF ARC PROFILES 


IF (FLOAT( INT(N/50)). EQ. FLOAT(N)/50. ) THEN 
WRITE(11) N,TIME, INT(YARC)-3,(TC(C(III, INT( YARC) -3,1)),III=1,27) 


ENDIF 


2. Source Listing of Subroutine FINLF Modifications 


the FIN subroutine. 


SUBROUTINE FINLF(B,C,COUT,LOF) 


‘ The top section of the subroutine is shown, the rest is identical to 


THIS SUBROUTINE HAS BEEN MODIFIED TO ALLOW INSERTION OF A LACK 

OF FUSION ZONE. THE ZONE MAY BE OF ANY SIZE BUT CANNOT INCLUDE 
POINTS WITH I<3 OR I>25 OR J<3 OR J>25 OR K>6 SINCE THESE BOUNDARY 
POINTS WOULD REQUIKE EXTENSIVE MODIFICATION OF THE SUBKOUTING: 

THE CALLING PROGRAM, WELDLF, POSITIONS THE LOF ZONE ON THE FINE 
GRID. JL SPECIFIES THE LENGTH ©F THE ZONE. I AND 12 SPECI a 
WIDTH OF THE ZONE. AND Kl AND K2 SPECIFY THE DEPTH OF THE ZONE. 


DIMENSION BC 27,27 510) 0@2 7227.8), COUI@27 27-6). OCs) ce hea 


DIMENSTON CK(27,27,8) 


COMMON NBSNC  YARC, SUM. TINE ao 


ty 


BOMMON / INPUT/ 11,12 ,K1,K2,90 
GK(T1,T2)=FK(T1)*(T1-T2) /(FKC(T1)+FK(T2)) 
DATA FKB/26. 5/ 


C CALCULATE THERMAL CONDUCTIVITY AT EACH NODE 
BO 10 T=1,27 
DO 10 J=1,27 
DO 10 K=1,8 
eM@ie I,K )=PRCGC1, Jk) ) 
10 CONTINUE 


C SET LACK OF FUSION ZONE TO ZERO CONDUCTIVITY 


DO 11 J=LOF-JL+1,LOF 
DO 11 I=11,I2 
DO 11 K=K1,K2 
11 CK(I,J,K)=0. 


* A. THE INTERIOR NODES 


DO 1 I=2,26 
DO 1 J=2,26 
DO 1 K=2,7 
COUT(I,J,K)=FK(C(1,J,K))*FO(3)*(HK(I+1,J,K,1,J,K,C,CK)+ 
* HK(I-1,J,K,1,J,K,C,CK)+HK(1I,J+1,K,1,J,K,C,CK)+HK(I,J-1,K, 
cae) CCK) +HK( 1,J3,K+1,1,3,k,C,CK)+HK(1,J,K-1,1,3,K,C,CK)) 
1 CONTINUE 


* 8B. THE TOP SURFACE NODES 


TINF4=TINF**4 
Be 2 1=2.26 
Den2 J=2.26 
COUT(I,J,1)=FK(C(I,J,1))*FO(3)*(HK(I+1,J,1,1,J,1,€,CK)+ 

meee Hk( I-1,3,1,1,J,1,C,CK)+HK(1,J+1,1,1,J,1,C,CK)+ 

fe HK(I,J-1,1,1,J,1,€,CK)+2. *HK(1,J,2,1,J,1,C,CK))+ 

* BI(3)*(TINF4-C(I,J,1)**4)+BI(4)*(TINF-C(I,J,1)) 
2 CONTINUE 


* The function HK, which performs save function as GK but uses the 
* precalculated thermal conductivities. 


REAL FUNCTION HK(I,J,K,L,M,N,C,CK) 
DIMENSION CK(27,27,8),C(27,27,8) 
A=CK(1I,J,K)+CK(L,M,N) 
IF (A.EQ. 0.) THEN 
HK=0. 
ELSE 
HK=CK(1,J,K)*(C(1,J,K)-C(L,M,N))/A 
ENDIF 
END 


100 FORMAT( 


3. Source Listing of Program TABLE 


PROGRAM TABLE 

REAL T020.27) .GbinCzo) 

INTEGER N(20),J(20) 

OPEN(1,FILE='COMP' ,FORM='UNFORMATTED' ) 
OPEN(2,FILE='PRINT ,STATUS='NEW' ) 

DO 1 K=1,20 


1 READ(1) NCK)o 2IMp GO rie), Clk) i) 


WRITE(2,100) (N(K),TIME(K),J(K) ,K=1,20) 
WRITE(2,200) (CINT((T(K,1)-300)) ,K=1,20) ,1=1,27) 
N TIME. ' YARC=1' (200 /14,88. 2 Js 


200 FORMAT( 2014) 


10 

20 

30 

40 

50 

60 

70 

80 

90 

100 
105 
IL 10 
120 
130 
140 
150 
160 
170 
180 
10, 
200 
210 
220 
230 
240 
Za0 
260 
270 
280 
290 
300 
310 
320 
330 
340 
350 
360 
370 
380 
3710 


END 
4. Source Listing of Program FLAW 


This program is used to process the data from running the program 
WELDLF. After running WELDLF, the program TABLE is run to produce 
the history of temperature profiles behind the arc for different 
geometries of flaws. This data is manually processed to determine 
the max temperature change from quasi-steadystate and the max 
percentage change from quasi-steadystate. This information, plus the 
flaw geometry and arc relative power is placed in a data file which 
can be read by this program in the order of: [1, 12, JL, Kiger 
Delta T max, Percent Delta T max, Relative Power. The program then 


‘performs fits the data to a relationship between the max percent 
‘change in temperature from quasi-steadystate and the location of 
the flaw. That relationship is: 

T - Tass CPA 

Tqss R*Z Zz 

where: 

X = (11+12)/2-14 

A = (12-11)*JL 

P = Relative Power 

Z = K1-1 

R = SQRT(X 2+ Z 2) 

C = Constant of the relationship 

CLS 

PRINT | This Program determines the Relationship between’ 
PRINT " Flaws and the Resultant Temperature Change." 
PRINT 

' Input the Flaw parameters 

INPUT "HOW MANY DATA POINTS -ARE IN THE FILE";N 

DIM X1(0N),X2(N),DYCN),Z10N) ,Z2(N), X(N), DXCN) GAG) ECND (DIC ia 


DIM P(N),Z(N),CT(N) 

INPUT "What is the name of the file’; FILES 

OPEN "I'' ,#1,FILES 

FOR I=1 TO N 

INPUT #1,X1(1),X2(1),DY(1), Z1@) 522 eae ene Cl) 
NEXT I 

' Echo the input 
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400 PRINT 

410 PRINT "No. ee 'f ee ENG ee Zale Oe, 

420 FOR I=1 TON 

PeeR INT J: ' :X1(1),X2C1),DY(1),2101),22C1) 

440 NEXT I 

450 ' Calculate the individual coefficient for each flaw. 
BeemPRINT : PRINT "No. ","AREA","R", "COR T" 

470 FOR I=1 TON 

480 DX(1)=ABS(X2(1)-X1(1)) +1 

490 X(I)=(X1(1)+X2(1))/2 - 14 

Byeez( 1)=Z1(1)-1 

510 RCI)=SQR(X(I1) 2+Z(1) 2) 

Boma 1 )=DX(1)*DY(1) 

530 CTC I)=P(1)*ACI)/CRCI)*2(1) 2) 

BogmPRINT 1,A(1),RC(1),CT(I) 

550 NEXT I 

560 INPUT "Strike Enter to continue"; KS 

570 ' Calculate the mean coefficient and display individual fits. 
580 PRINT 

590 PRINT "No.", "TEMP" ," PERCENT" ,"TEMP COR",'PERCENT COR" 
600 FOR I=1 TON 

GMOmERINT I ,DI(1),PICL),DT(1)/CTC(1),PT(1)/CTCI) 

eee I—PT(1)/CTC1) 

630 SUM=SUM+XI 

640 S2=XI*XI+S2 

650 NEXT I 

660 S=SQR((S2-SUM*SUM/N)/(N-1)) 

670 PRINT : PRINT “AVERAGE CORRELATION FIT = ";SUM/N,"STD DEV = ";S 
680 END 


5. Listing of a typical exec for repeated program execution 


& TRACE 

COPYFILE BASIC * B FILE = A 
&STACK 11,14,4,4,3 

EXEC RUN WELDLF 

EXEC TDISK 15 DIS 

EXEC DISSPLA PLOT1 

EXEC SHERPA TEMPOUTP SHGRAPH T 
EXEC DISSPLA PLOT2 

EXEC SHERPA TEMPOUTP SHGRAPH T 
EXEC DISSPLA PLOT3 

EXEC SHERPA TEMPOUTP SHGRAPH T 
EXEC PRINT FILE VERIF A 
COPYFILE FILE COMP A = COMP9 B 
EXEC DISSPLA PLOTC 

EXEC SHERPA TEMPOUTP SHGRAPH T 
EXEC RUN TABLE 

EXEC PRINT FILE PRINT A 

ERASE FILE * A 

COPYFILE BASIC * B FILE = A 
emack 11,14,4.4,2 

EXEC RUN WELDLF 

EXEC DISSPLA PLOT1 

EXEC SHERPA TEMPOUTP SHGRAPH T 
EXEC DISSPLA PLOT2 

EXEC SHERPA TEMPOUTP SHGRAPH T 


EXEC DisskLA PLOTS 

EXEC SHERPA TEMPOUTP SHERAriIat 
EXEC PROMS ILE VERTE 4 
COPYFILE FILE COMP A = GOMPIOrE 
EXEC DISSPLA PLOTC 

EXEC SHERPA TEMPOUTP SHGRAPH T 
EXEC RUN TABLE 

EAEC ERING PILES PRINT of 

ERASE a Die * 74 

FE. COOLING RATE PROGRAMS 

To measure the cooling rates during startup, steady state and shutdown the program 
WELD was modified and called WELDC. The modified subroutines are called FINC, 
MEDC and CORC. An additional subroutine RATES 1s added for calculating the 
cooling rates. Since the temperatures of interest are further from the arc, the Fine and 
Medium zones were extend to allow measuring these temperature. This modification 
approximately doubles the time to run the program. The other major modification was 
the wav in which the heat is added to the material and how the molten region 1s modeled. 
The heat is only added to the surface nodes. To simulate the effect of arc pressure, a 
directional thermal conductivity is used, where the thermal conductivity in the Z direc- 
tion 1s ten times that in the horizontal directions. This later effect only applies to those 
nodes which are above the melting temperature of the metal. This was implemented by 
an additional version of difference function GK, called HK, which is only used for dif- 
ferencing in the Z direction. 

The RATES subroutine is called once every time step. It measures three separate 
cooling rates for the surface nodes in the heat affected zones. These cooling rates are 
determined by measuring the time it takes the temperature at a given absolute position 
to cool through a specified temperature range. This is performed in a three step process 
over all of the nodes currently in these fine and medium areas behind the arc. It first 
checks to see if the node temperature is now greater than the upper limit of the tem- 
perature band. If it 1s, this 1s noted for future reference with the logical array CT. The 
second check is to see if a node temperature is below the upper limit and previously 
above the upper band and vet never below it before. If this is the first time it has re- 
turned to below the upper band. this fact is noted for future reference in the logical array 
CT2. In addition, the time and temperature for this node are saved in the array CR. 


The third check is to see if a node temperature is below the lower limit of the temper- 


ature band having been previously above it. If it has returned to below the lower limit 
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for the first time then this is noted by setting the array CT for that node to false and the 
time and temperature for this node is again saved in the array CR. 

This procedure is done in both the fine and medium zones. Due to the difference in 
sizes, Only one out of nine fine nodal points will have a corresponding medium node. 
Thus if a nodal point temperature is not below the lower limit before the fine grid passes 
on, the program will not always report a cooling rate at the node. The subroutine 
RATES passes the arrays back to the main program WELDC which saves these files so 
that the program can be restarted. In addition, the CR arrav 1s separately outputted so 
that it can be processed to determine the cooling rates. 

The program that determines the cooling rates is called RATEOUT. It reads the 
data file produced by WELDC and converts the data to cooling rates. It also provides 
the Y and R location of the node from the arc position at the mid-time of the cooling 
interval, as Well as the absolute location of the node. If no cooling rate was determined 
for a location it will then provide some indication of the status of that location. 

Due to the extensive modifications required to implement these changes, the new 
version of the source code is included in full. Extreme care was required to change all 
of the indices, it Was almost as much work as writing the original code. 

1. Source Listing of Program WELDC 


PROGRAM WELDC 


sedvtevoovvlesevevedesevedsdouevedesesesedesevescdesvacslevesevedescoevedevedvscvestededeved deskcaevedaedevedekesecsvedvedcvesevevesevevevevevevese 


re ve 
* DATE: 16 AUGUST 1988 WRITTEN BY: ROBERT ULE ¥ 
* DATE: 17 AUGUST 1988 MODE TEDe By: SROBERE ULE ¥ 
* USES DIRECTIONAL THERMAL CONDUCTIVITY IN THE WELD POOL AND % 
* LONGER ARRAYS IN THE J DIRECTION TO STUDY COOLING RATES. ¥ 


ve THE ARC IS ADDED ONLY TO THE SURFACE NODES. ve 
Sel de else geese ote ae ae \ er) ed Vohra ae ob ke ek oe) he) de ee oe ee ee er) oe) te) ek Se) Ne hed) le he i oh) ki) lh od li 
DIMENSION A(21,36),AOUT( 21,36) ,AIN( 21,36) ,ASUM( 21, 36), 
checyenscrlO) .b000( 27 .36,10),BINC 27 , 36,10), BSUM( 27 , 36,10), 
*C(27,45,8) ,COUT(27,45,8),CIN(27,45,8) ,CSUM(27,45,8), 
eros he hl¢4),CR(200,14,3,4) ,CTC200, 14,3) 
INTEGER NDIV,I,J,K,M,BSTEP,CSTEP 
CHARACTER ANS 
POGLeAINY ES .CT,C1T2( 200, 14,3) 
COMMON NB,NC,YARC,SUM,TINF,BI,FO 


DATA A,B/10476**300. 0/,C/9720%1. 14237E8/,ASUM, BSUM/10476*0. 0/ 
% AOUT , BOUT/104767*0. 0/,CT/8400*. FALSE. /,CT2/8400* 
*, FALSE. / ,CSUM/9720**0. /,CR/33600*0. / 
C SET TIME TO RUN THE PROBLEM AND FIXED CONDITIONS 


FINI = 20. 
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NDIV=FINI*100 
XVEL=. 0 
DELT=. 01 
TINF=300. 0 


IF THIS IS A RESTART OF A PREVIOUS PROBLEM GOTO 100 
GOTO 100 
SET INITIAL VALUES FOR STARTING A PROBLEM 


YVEL=4. 
VOLT=30. 
AMP=100. 
EPP —1 65 


OPEN Ibe OUTPUISPITES 


OPEN(1,FILE='SURF' ,STATUS=' NEW’ , FORM=' UNFORMATTED' ) 
OPEN(2,FILE='FINAL’ ,STATUS='NEW , FORM='UNFORMATTED' ) 
OPEN( 3 ,FILE='CUT' , STATUS='NEW' , FORM='UNFORMATTED' ) 

OPEN(4,FILE='RATE ,STATUS='NEW’ ,FORM='UNFORMATTED' ) 


THE INITIAL CONDITIONS 


TIME=0. 
STEP=3. 
OUT=. 499 
N=0 
BSTEP=0 
CSTEP=0 
NB=3 
NC=10 


WELD PARAMETERS 


QDOT=EFF*VOLT*AMP 
Q=QDOT*DELT*1. E9 


REENTRY FOR RESTART 
200 CONTINUE 
BOUNDARY CONDITION COEFFICIENTS 


FO( 1)=DELT*. 1636 
FO(2)=DELT*1. 4722 
FO(3)=DELT*1000000. 
BIG) =s00132 

BIg 2) =s00T132 
BI(3)=DELT*. 00009299 
BI(4)=DELT*50000. 


CALCULATE THE RUNGE-KUTTA APPROXIMATION 


DIS=TIME*YVEL 
YARC=YVEL* ( TIME+DELT/2)+91-NB*9 -NC*3 
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DO 10 M=1,NDIV 
TIME=TIME+DELT 
DIS=TIME*YVEL 
N=N+1 


C POSITION HEAT SOURCE AND CALCULATE VOLUME WEIGHTING FACTOR  SUM' 


YARC=YVEL**( TIME+DELT /2)+9 1-NB*9 -NC*3 
XARC=14. 0 
SUM=0. 
DO 1 J=28,41 
DO 1 I=9,19 
SUM=SUM+EXP( -. 10625%( ( J-YARC)**2+( I-XARC)**2))*. 5 
DO 1 K=1,3 
SUM=SUM+EXP( -. 10625**( ( J-YARC)**2+( I-XARC) **2+K%**2 ) ) 
1 CONTINUE 
SUM=SUM/Q 


C ADD THE HEAT FROM THE ARC 


Hong J=25 901 
HOmze=9), 19 
C(I,J,1)=C(1,J,1)+EXP(-. 10625*( (J-YARC)**2+( I-XARC)**2)) /SUM 
*, 5 


DOe2aK=2ne 
C(CI,J,K)=C(1,J,K)+EXP( -. 10625*( ( J-YARC)**2+( I -XARC)***2 
= +(K-1)**2))/SUM 
2 CONTINUE 


C CONVERT ENTHALPY TO TEMPERATURE FOR THE FINE PLATE 


BO 60 I=1,27 
DO 60 J=1,45 
DO 60 K=1,8 
CANGIek)=1( CC 1,.J4K)) 
60 CONTINUE 


C CALCULATE COOLING RATE DATA 
CALL RATES(B,CIN,CR,CT,CT2,NB,NC,TIME) 
C FIRST RUNGE-KUTTA ITERATION 
CALL CORC(A,B,ASUM) 
CALL MEDC(A,B,CIN, BSUM) 
CALL FINC(B,CIN,CSUM) 
C FIND T+K1/2 
DO 11 I=1,21 
DO 11 J=1,36 
11 AIN(I,J)=A(I,J)+ASUM(I,J)/2. 
DO 13 I=1,27 


DO 13 J=1,36 
DO 13 K=1,10 
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13. BIN(I,J,K)=B(1,J,K)+BSUM(1I,J,K)/2. 
DO 12 I=1,27 
peri? J=1.45 
BO) 2 hl 
CIN(I,J,K)=T(C(I,J,K)+CSUM(1,J,K)/2. ) 
12 CONTINUE 


SECOND ITERATION 


CALL CORC(AIN, BIN, AOUT) 
CALL MEDC(AIN, BIN,CIN, BOUT) 
CALL FINC(BIN,CIN,COUT) 


FIND T+K2/2 


DO. 27, Ten 
DO 21 J=1,36 
ASUM( I, J)=ASUM(I, J) +2. *AOUT(I, J) 
21 AIN(I,J)=A(1,J)+AOUT(I,J)/2. 
We, BE WH Dy 
DO 23 J=1,36 
DO 23 K=1,10 
BSUM(1,J,K)=BSUM(1I,J,K)+2. *BOUT(I,J,K) 
23. BIN(1I,J,K)=BC(I,J,K)+BOUT(I,J,K)/2. 
0) 2D Wei ay 
DO 22 J=1,45 
DO 22 K=1,8 
CSUM(I,J,K)=CSUM(I,J,K)+2. *COUT(I,J,K) 
CIN(I,J,K)=T(C(1,J,K)+COUT(1I,J,K)/2. ) 
22 CONTINUE 


THIRD ITERATION 
CALL CORC(AIN, BIN, AOUT) 


CALL MEDC(AIN, BIN, CIN, BOUT) 
CALL FINC(BIN,CIN,COUT) 


FIND T+K3 
DO 31 I=1,21 
DO 31 J=1,36 


ASUM( I, J)=ASUM(I, J)+2. *AOUT(I, J) 
31 AIN(I,J)=A(I,J)+A0OUT(I, J) 
Dors3 I=127 
DO 33 J=1,36 
DO 33 K=1,10 
BSUM(I,J,K)=BSUM(1I,J,K)+2. *BOUT(1I,J,K) 
33. BIN(I,J,K)=B(1,J,K)+BOUT(I,J,K) 
DO 32 I=1,27 
DO 32 J=1,45 
DO 32 K=1,8 
CSUM(I,J,K)=CSUM(I,J,K)+2. *COUT(I,J,K) 
INCI. J K)=I( CCl ink) +COUm@IeeIa an 
32 CONTINUE 


FORTH ITERATION 
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CALL CORC(AIN, BIN, AOUT) 
CALL MEDC(AIN,BIN,CIN, BOUT) 
CALL FINC(BIN,CIN, COUT) 


PERFORM THE RUNGE-KUTTA SUM 


DO 40 I=1,21 
DO 40 J=1,36 
40 A(1,J)=A(1,J)+(ASUM(I,J)+AOUT(I,J))/6 
DO 43 I=1,27 
DO 43 J=1,36 
DO 43 K=1,10 
43 B(1,J,K)=B(1,J,K)+(BSUM(I,J,K)+BOUT(I,J,K))/6 
DO 41 I=1,27 
DO 41 J=1,45 
DO) il eal os 
mee GC 1,J,K)=C(1,J,K)+(CSUM(1 ,J,K)+COUT(1,J,K))/6 


ser GRID IF REQUIRED 


BP yODIs.GE.STEP) THEN 
STEP=STEP+3. 
BSTEP=BSTEP+1 


MOVE FINE GRID, FIRST AVERAGE FINE ENTHALP AND ASSIGN THE EQUIVALENT 
TEMPERATURE TO THE MEDIUM GRID 


DO 51 IB=1,9 
HB1=0. 
HB2=0. 
HB3=0. 
DO 50 IC=IB*3-2, IB*3 
DO 50 JC=1,3 
HB1=C(IC,JC,1)+2.*C(1IC,JC,2)+HB1 
DO 50 KC=3,5 
HB2=HB2+C(IC,JC,KC) 
HB3=HB3+C(IC,JC,KC+3) 
50 CONTINUE 
B( IB+9 ,NC,1)=T(HB1/27) 
B( IB+9 ,NC, 2)=T(HB2/27) 
B( IB+9 ,NC, 3)=T(HB3/27) 
51 CONTINUE 


SHIFT FINE GRID 
DO 52 J=1,42 
JI=I+3 
DO 52 K=1,8 
Mem52 1=1,27 
oo Ge ek)=C(1,JJ,K) 
ASSIGN MED TEMPS TO FINE GRID 


DOws> 7 1=10,18 
TB1=B( I ,NC+15,1) 
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TB2=B(1,NC+15,2) 
TB3=B(1,NC+15,3) 
DO 53 IC=(1-9)*3-2,(1I-9)*3 
DO 53 JC=43,45 
C(IC,JC,1)=H(TB1) 
C(IC,JC,2)=H(TB1) 
DO 53 KC=3,5 
C( IC, JC,KC)=H(TB2) 
C(IC,JC,KC+3)=H(TB3) 
53 CONTINUE 
NC=NC+1 


CHECK TO SEE IF TIME TO SHIFT MEDIUM Grip 
IP (BSTEP. Bato acy 
BSTEP=0 


THEN AVERAGE MEDIUM GRID AND ASSIGN TO COARSE GRID 


DO 54 IA=1,9 
WME. 
DO 55 IB=IA*3-2,IA*3 
DO 55 JB=1,3 
TA1=TA1+(B(IB,JB,1)+B(1IB,JB,10))/2 
DO 55 KB=2,9 
TA1=TA1+B(1B,JB,KB) 
55 CONTINUE 
A( IA+6 ,NB)=TA1/81 
54 CONTINUE 


SHIFT MED GRID 


DO 56 J=1,33 
JI=I+3 
DO 56 I=1,27 
DO 56 K=1,10 
56 B(I,J,K)=B(1,JJ,K) 


ASSIGN COARSE TEMPS TO MEDIUM GRID 


DO 57 1=7.15 
TA1=A(I ,NB+12) 
DO 57 IB=(1I-6)*3-2,(1-6)*3 
DO 57 JB=34, 36 
DO 57 KB=1,10 
B(IB,JB,KB)=TA1 
57 CONTINUE 
NB=NB+1 
NC=10 
ENDIF 
ENDIF 


OUTPUT RUNNING SOLUTION EVERY .5 SECONDS 


LE eG eivbeore OUI). tin 
OUT=O0UT+. 5 


Paley ( CAMGGCT Jil), allan en ees 5) 
AVERAGE FINE ENTHALPIES AND ASSIGN EQUIVALENT TEMPERATURE TO THE 
MEDIUM GRID 

DO 61 J=0,14 

DO 61 IB=1,9 
HB1=0. 
DO 66 IC=IB*3-2,1B*3 
DO 66 JC=J*341,J*3+3 
66 HB1=C(IC,JC,1)+2*C(1C,JC,2)+HB1 
61 B(9+IB ,NC+J,1)=T(HB1/27. ) 
WRITE(1)((B(I,J,1) ,I=1,27) ,J=1, 36) 
WRITE(3)((TCC(1, INT( YARC) ,K)) ,I=1,27),K=1,8) 
ENDIF 
10 CONTINUE 


OUTPUT FINAL RESULTS 


WRITE(2) TIME 
WRITE(2)((( T(CCI,J,K)),1=1,27),J=1,45) ,K=1,8) 
AVERAGE FINE ENTHALPIES AND ASSIGN EQUIVALENT TEMPERATURE TO THE 
MEDIUM GRID 
DO 62 J=0,14 
DO 62 IB=1,9 
HB1=0. 
HB2=0. 
HB3=0. 
DO 63 IC=IB*3-2,1B*3 
DO 63 JC=14+J*3 , 3433 
HB1=C(IC,JC,1)+2.*C(IC,JC,2)+HB1 
DO 63 KC=3,5 
HB2=HB2+C(IC,JC,KC) 
HB3=HB3+C( IC, JC, KC+3) 
63. CONTINUE 
B( IB+9 ,NC+J, 1)=T(HB1/27. ) 
BC IB+9 ,NC+J, 2)=T(HB2/27. ) 
BC IB+9 ,NC+J ,3)=T(HB3/27. ) 
62 CONTINUE 
WRITE(2)(((B(1,J,K) ,I=1, 27) ,J=1,36) ,K=1, 10) 
AVERAGE MEDIUM TEMPERATURES AND ASSIGN TO COARSE GRID 
DO 64 J=0,11 
DO 64 IA=1,9 
TA1=0. 
DO 65 IB=IA*3-2,IA*3 
DO 65 JB=14+J*3, 3433 
TA1=TA1+(B(IB,JB,1)+B(IB,JB,10))/2 
DO 65 KB=2,9 
TA1=TA1+B( IB, JB,KB) 
65 CONTINUE 
A( IA+6 ,NB+J)=TA1/81 
64 CONTINUE 
WRITE(2)((ACI,J),1=1,21),J=1, 36) 


ALLOW OPTION OF RESTARTING THE PROBLEM 


INQUIRE(9 ,OPENED=YES ) 
TEC CES ea HEN 


Is2 


REWIND(9 ) 
ELSE 
OPEN(9,FILE='RESTAR' ,STATUS='NEW' , FORM='UNFORMATTED' ) 
ENDIF 
WRITE(9)OUT,Q,NB,NC,YVEL,STEP,CSTEP,BSTEP,N,DELT 
WRITE(9)A 
WRITE(9)B 
WRITE(9)C 
WRITE( 4) GR C202 
STOP 


CG PROCEDURE FOR RESTARTING A PREVIOUS PROBLEM 


100 OPEN(2,FILE='FINAL' , STATUS='OLD' , FORM='UNFORMATTED' ) 
OPEN(4,FILE='RATE' ,STATUS='OLD' , FORM='UNFORMATTED' ) 


C THE INITIAL CONDITIONS 


READ( 2) TIME 

REWIND( 2) 

IF (TIME. LT. (.499)) THEN 
OPEN(1,FILE='SURF' , STATUS=' NEW' , FORM='UNFORMATTED' ) 
OPEN(3,FILE='CUT' ,STATUS='NEW' , FORM='UNFORMATTED' ) 

ELSE 
OPEN(1,FILE='SURF' ,STATUS='OLD' , FORM='UNFORMATTED' ) 
OPEN(3,FILE='CUT' ,STATUS='OLD' ,FORM='UNFORMATTED' ) 

ENDIF 

OPEN(9,FILE='RESTAR’ ,STATUS='OLD' , FORM='UNFORMATTED' ) 

READ(9)OUT,Q,NB,NC,YVEL,STEP,CSTEP,BSTEP,N,DELT 


C POSITION THE RUNNING OUTPUT FILES TO THE END OF THE FILES 


DO 102 L=1,N/50 
READ(1)((C(1,J,1),I=1,27),J=1,45) 
READGIIGCCE( aint) 1=1.27) I=) co) 
READ(3)((C(1,1,K),I=1,27) ,K=1,8) 

102 CONTINUE 

READ(G)UGR- CT. CT2 enh 
REWIND(4) 


C INPUT THE TEMPERATURE/ENTHALPY MATRICES 


READ(9)A 
READ(9 )B 
READ(9)C 
GOTO 200 
END 


2. Source Listing of Subroutine FINC 
SUBROUTINE FINC(B,C,COUT) 
C THIS SUBROUTINE IS ASSOCIATED WITH THE PROGRAM WELDC, WHICH USES A 


C DIRECTIONAL THERMAL CONDUCTIVITY IN THE MOLTEN POOL. IT HAS EXTENDED 
C ZONES IN THE J DIRECTION TO ALLOW BETTER RESOLUTION OF COOLING RATES. 


DIMENSION B(27,36;10) ,C(27,45,8) ,COUI( 27 aaa Fe lca 
COMMON NB NC, YARC,5UM TINE Bis bo 
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en Gee, 2) =FKCT1 (11-712) /( FK(1T1)+FKC22) ) 
DATA FKB/26. 5 
C THE INTERNAL NODES 


DO 1 I=2,26 
DO 1 J=2,44 
DO 1 K=2,7 
COUT(1,J,K)=FK(C(I, 1J,K) *FOC3)*(GK(C(I+1, 3 pone oo 
* fe en) ett .K) )+GK(C(I,J#1,K) ,CC1,J,K))+GK(C(I,J-1,X), 
oemcc me oh) )tHRCC(L,J,K+1) ,CC1, J ,K) )+HK(C(I,J hee CGigwN) )) 
1 CONTINUE 


C THE TOP SURFACE 


TINF4=TINF* 4 
DO 2 I=2,26 
DO 2 J=2,44 
COUT(1,J,1)=FK(C(I,J,1))*FO(3)*(GK(CC(I+1,J,1),C(1,J,1))+ 
Beer eeei- sg e1) OCT. J, 1))4+GK(C(1,J+1.1).C(1,3,1))+ 
ba Meese) 0-1 1). C(1,J, 1) )4+2. =HK(C(1,J,2),C(1,J,1)))+ 
Be 3) (PINES -C( 1) J, 1)**4)+B1(4)4(TINF-C(I,J,1)) 
2 CONTINUE 


GC THE BOTTOM SURFACE (MEDIUM INTERFACE) 


Demer i=? 26 
IB=(I-. 5)/3+10 
BO 3 J=2,44 
JB=(J-.5)/3+NC 
POU 6 =F 0( 3) =( 7 K(CC1,J, 8) )*CGR(CCI4+1,J,8),C(1,J,8))+ 
Emer CCC I-13 .8),CCl,J,8))+GK(C(1,J+1,8),C(1,J,8))+GK(C(1,J-1,8), 
* © C(1,J,8))+HK(C(1,J,7),CC1,J,8)))+FKB*(BCIB,JB,4)-C(1,J,8))) 
3 CONTINUE 


GC THE ENDS (MEDIUM INTERFACE) 


JB=NC-1 
JMAX=NC+15 
DO 4 I=2,27 
IB=(I-.5)/3+10 
DO 4 K=2,7 
KB=2 
IF(K.EQ.2) KB=1 
IF(K. GE. 6) KB=3 
GOUT(I,1,K)=FO(3)*(FK(C(1I,1,K))*(GK(C(1+1,1,K),C(1,1,K))+ 
ec een et onc 113k) )+CKCGC 1, 23K), €C1, 15K) )+GK(CC1,1,K+1), 
*% (C(1,1,K))+GK(C(1I,1,K-1),CC1I,1,K)))+FKB*(BCIB,JB, KB) - 
eGipete),)) 
COUT(I,45 ,K)=FO(3)*(FK(C(1,45,K))*(GK(C(1+1,45,K),C(1,45,K))+ 
# GK(G(1-1,45,K),C(1,45,K) )+GK(C(1,44,K),CC(1,45,K))+ 
*  HK(C(1I,45,K+1),C(1,45,K) )+HK(C(1,45,K-1),C(1,45,K)))+ 
* -FKB*(BC IB, JMAX,KB)-C(1,45,K))) 
4 CONTINUE 


st. 
“ 


[5° 


C THE 


5 


C THE 


6 


C THE 


SIDES (MEDIUM INTERFACE) 


DO 5 J=2,44 
JB=(J-.5)/3+NC 
DO 5 K=2,7 
KB=2 
IF(K. EQ. 2) KB=1 
IF(K. GE. 6) KB=3 
COUT( 1,J,K)=FO( 3)*(FK(C(1,J,K))*(GK( GC 2.01, K)) Gq eee nies 
GK(C(1,J+1,K),G(1,4,K))+GR0GC1,J-1k)), cc aone 
HK(C(1,J,K+1),C(1,J,K))+HK(C(1,J,K-1),C(1,J,K)))+ 
FKB*(B(9,JB,KB)-C(1,J,K))) 
COUT( 27,J,K)=FO(3)*(FK(C(27,J,K))*(GK(CC26,J.K) C027 sie hones 
GK(G(27,J+1,K) Gl 27, KGR( C27, J- 1 emcee Toa 
HK(C(27,J,K+1),C(27,J,K) )+HK(C(27,J,K-1),C( 27 aeons 
FKB*(B(19,JB,KB)-C(27,J,K))) 
CONTINUE 


END EDGES (MEDIUM INTERFACE AND SURFACE EFFECTS ON TOP) 


KB=1 
JB=NC-1 
JMAX=NC+15 
DO 6 I=2,26 
IB=(I-.5)/3+10 
COUT( 1,1, 1)=FO(3)*( FK(C(1,1,1))*(GK(C(1+1,1,1),C(1,1,1))+ 
GK(C(1-1,1,1),CC1,1, 1))+GKGGer 2) nec a een 
2. *GKCO(1 21.2) .0( 1.111) ))aPRE* (BC IE Re KepeG( Ele er 
BI(3)*(TINF4-C(1,1,1)**4)+BI(4)*(TINF-C(I,1,1)) 
COUT(1,45,1)=FO(3)*(FK(C(1,45,1))*(GK(C(I+1,45,1),C(1,45,1))+ 
GK(C(1-1,45,1).C(1.45,1))+GK(GC 1. 24. 1 nO Gi Seine: 
2. *GK(C(1,45,2),C(1,45,1)))+FKB*(B( IB, JMAX ,KB)-C(1I,45,1)))+ 
BI(3)*(TINF4-C(1,45,1)**4)+B1I(4)*( TINF-C(1,45,1)) 
COUT(1I,1,8)=FO(3)*(FK(C(1,1,8))*(GK(C(1+1,1,8),C(1,1,8))+ 
GKCCCI-1 5,158) 56Ci. 155) ener 25> Olen ens 
GK(C(I,1,7),C(1,1,8)))+FKB*(B(IB,JB,3)+B(IB,NC,4) 
22 C1, 1.8) 
COUT(1,45,8)=FO(3)*(FK(C(1I,45,8))*(GK(C(I+1,45,8) ,C(1,45,8))+ 
GK(G(1-1,4558) .C( 1.45.6) 4CRCGGL 64) GG Cnecseepe 
GK(C(1,45,7) ,C(1,45,8)) )+FKB*(B( IB, JMAX,3)+B( IB,NC+14,4)- 
2.*C(1,45,8))) 
CONTINUE 


SIDE EDGES (MEDIUM INTERFACE AND SURFACE EFFECTS ON TOP) 


DO 7 J=2,44 

JB=( I-25 )/ 3+NG 

COUT( 1,J,1)=FO( 3)*(FK(C(1,J,1))*(GK(C(1,J+1,1),C(1,J,1))+ 
GKCCC 1, J-15 1) iGGi5 J, 1) ) +GCKCC C2 ee) CC ene 
2.*GK(C(1,J,2),C(1,J,1)))+FKB*(B(9,JB,1)-C(1,J,1)))+ 
BI(3)*(TINF4-C(1,J,1)**4)+B1(4)*(TINF-C(1,J,1)) 

COUT( 27 ,J,1)=FO(3)*(FK(C(27,J,1))*(GK(C( 27 ,Jt1,1),C(27,J,1))+ 
GK(C(27,J-1,1) ,C(27,J,1))+GK(C(26,J,1),C(27,J,1))+ 
2.*GK(C(27,J,2),C(27,J,1)))+FKB*(B(19,JB,1)-C(27,J,1)))+ 
BI(3)*(TINF4-C(27,J,1)**4)+B1I(4)*(TINF-C(27,J,1)) 
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COUT(1,J,8)=FO(3)*(FK(C(1,J,8))*(GK(C(1,J+1,8),C(1,J,8))+ 

= GRCEGlC I-18), CCl ne )) +GNCGC24 6), 60 158) + 

*  GK(C(1,J,7),C(1,J,8)) )+FKB*(B(9,JB,3)+B(10,JB,4) 

*  -2.%C(1,3,8))) 
COUT(27,J,8)=FO(3)*(FK(C(27,J3,8))*(GK(C(27,J+1,8) ,C(27,J,8))+ 

Peo Ceara ons G2 ane ont OC2o,dsc) ,6(27,J4,6) )+ 

* GK(C(27,J,7),€(27,J,8)))+FKB*(B(19,JB,3)+B(18,JB,4) 

* -=2,*C(27,3,8))) 

7 CONTINUE 


GC THE CORNER EDGES 


Doms eK? 7 
KB=2 
IF(K.EQ. 2) KB=1 
IF(K.GE.6) KB=3 
Counc 1h) =EOGs )(ERGCC dak) *CGK(C(2,1,),0¢1,1,k))- 
eo GC (dn creel hho ron C@l ll Kel) CCl ak) )o 
*  GK(C(1,1,K-1),C(1,1,K)))+FKB*(B( 10,NC-1,KB)+B(9,NC,KB)- 
pee 61,1 .K))) 
COUT(1,45 ,K)=FO(3)*(FK(C(1,45,K))*(GK(C(2,45,K),C(1,45,K))+ 
*  GK(C(1,44,K),C(1,45,K))+GK(C(1,45,K+1),C(1,45,K))+ 
*  GK(C(1,45,K-1),C(1,45 ,K)))+FKB*(B(10,NC+15 ,KB)+B(9,NC+14,KB)- 
ee 1 G5UK))) 
oume27 91K Y="OC 3)2*( FRC C( 27, 1,.K))*CGK(G( 26, 1K) .€(27,1,K))+ 
PC QGeie 7 2K) 027, VK) )tCK( (27 -1,K41) .6027,1,K) 
*  GK(C(27,1,K~-1),C(27,1,K)))+FKB*(B(18,NC-1,KB)+B( 19 ,NC,KB)- 
Same? (27. 1.K))) 
COUT( 27,45 ,K)=FO(3)*(FK(C(27,45,K))*(GK(C(27,45,K) ,C(27,45,K))+ 
pee RCGG27 44.8). C027 , 45. K) )+GK( C027, 45 .K41)),6027 .45.,K))+ 
* GK(C(27,45,K-1) ,C(27,45,K)) )+FKB*(B(19,NC+14, KB)+ 
*  B(18,NC+15,KB)~-2. *C(27,45,K))) 
8 CONTINUE 


C THE CORNERS (MEDIUM INTERFACE AND TOP SURFACE EFFECTS) 


Goumel. dl y=—O( SC FR C( 1.1. 1))*(GR(C(2,1.1) 9601 ,1,1))+ 

PCC O ll 2 wlan lad tL tGkh( CCl. 2) OC 1.151) ) J TFKB* 

* (B(10,NC-1,1)+B(9,NC,1)-2.*C(1,1,1)))+ 

* BI(3)*(TINF4-C(1,1,1)%"4)+BI(4)*(TINF-C(1,1,1)) 
G@owme 1,45, 1)=F0(3)*(FK(C(1, 45, 1))*(GK(C(2,45,1),C(1,45,1))+ 

peer KC OCW. 4201) .0( 1n45. 1) )+CKCECI.45 .2).0(1,45,1)) )+EKB* 

*  (B(10,NC+15,1)+B(9,NC+14,1)-2. *C(1,45,1)))+ 

* «BI 3)*(TINF4-C(1,45, 1)%"4)+B1I(4)*(TINF-C(1,45,1)) 
Geum 27 ,1,1)=FO(3)*(FRCC( 27, 1,1) )*CGK(C(26,1,1),C(27,1,1))+ 

GaGa 7 12 .1).G(27,1,1) )4GK0G( 27, 1.2),6( 27,1, 1) ))+FKB* 

*% (BC18,NC-1,1)+B(19,NC,1)-2.*C(27,1,1)))+ 

* BIC3)*(TINF4-C(27,1,1)**4)+BI(4)*(TINF-C(27,1,1)) 
COUT(27,45,1)=FO(3)*(FK(C(27,45,1))*(GK(C(26,45,1),C(27,45,1))+ 

poe Game? 7 4h. 11),,0( 27.45, 1) )+GK( (27,45, 2), 60 27,45.1)))+FKB* 

*  (B(18,NC+15,1)+B(19,NC+14,1)-2.*C(27,45,1)))+ 

*% BI(3)*(TINF4-C(27,45,1)*%4)+BI(4)*(TINF-C(27,45,1)) 
COUT(1,1,8)=FO(3)*( FK(C(1,1,8))*(GK(C(2,1,8),C(1,1,8))+ 

foe IGKGOGll 26) OC 1h 1.8) )4GKCGC 141-7). CC 1.1 eye KB* 

* (B(10,NC-1,3)+B(9,NC,3)+B(10,NC,4)-3. *C(1,1,8))) 
COUT(1,45,8)=FO(3)**(FK(C(1,45,8) )*(GK(C(2,45,8),C(1,45,8) )+ 
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*  GK(C(1,44,8),C(1,45,8))+GK(C(1,45.7)2601,45 60) ) +EKEs 

*  (B(10,NC+15,3)+B(9,NC+14,3)+B(10,NC+14,4)-3. *C(1,45,8))) 
COUT( 27,1,8)=FO(3)*CFK(C( 27, 1,8))*CGRCC( 26, 1,8), 0027.1 

*  GK(C(27,2,8),6(27,1,8))+GK(C(27,1,7).,c(27 incre: 

*  (B(18,NC-1,3)+B(19,NC,3)+B(18,NC,4)-3. *C(27,1,8))) 
COUT( 27,45 ,8)=FO(3)*(FK(C(27,45,8))*(GK(C(26,45,8) ,C(27,45,8))+ 

*% GK(C(27,44,8),C(27,45,8) )+GK(C(27 ,45,7),C(27,45,8) ))+EKB 

* (BC 18,NC+15,3)+B(19,NC+14,3)+B(18,NC+14,4)-3. *C(27,45,8))) 


RETURN 
END 


C FUNCTION FOR FINDING THE EFFECT OF THERMAL CONDUCTIVITY WHEN DIRECT- 
C DIRECTIONAL IN THE WELD POOL. 


FUNCTION HK(T1,T2) 
A=FK(T1) 
B=FK(T2) 
C=1. 
IF (A.EQ.35.) A=175. 
IF (B.EQ.35.) THEN 
B=175. 
c=5. 
ENDIF 
HK=A*(T1-T2)*C/( A+B) 
END 


3. Source Listing of Subroutine MEDC 


SUBROUTINE MEDC( A,B,C, BOUT) 
C THIS VERSION HAS BEEN MODIFIED TO HAVE LONGER ZONES IN THE J 
C DIRECTION TO ALLOW THE STUDY OF COOLING RATES. 


DIMENSION A(21,36),BOUT(27,36,10),B(27,36,10),BSUM(27,36,10), 
*C(27,45,8) ,FO(3),BI(4) 
COMMON NB,NC,YARC,SUM,TINF,BI,FO 


* A. THE INTERIOR NODES 


DO l=2 6 
DO) ged=2ees 
DO 777 K=2,4 
IF((1. GE. 9. AND. I. LE. 19. AND. J. GE. NC-1. AND. J. LE. NC+15) 
. AND. . NOT. (K. EQ. 4. AND. (I. EQ. 9. OR. I. EQ. 19. OR. J. EQ. NC-1. OR. J. EQ. 
* NC+15)). AND. . NOT. ((J. EQ. NC-1. OR. J. EQ. NC+15). AND. (1. EQ. 9. OR. I 
+ WECM como gn 
BOUT(1,J,K)=FO( 2)*(B(I+1,J,K)+B(I-1,J,K)+B(1I,J-1,K)+B(1I,J+1,K)+ 
* B(1,J,K+1)+B(1,J,K-1)-6. *B(1,J,K)) 
aa CONTINUE 
DO 7 K=5,9 
BOUT(1,J,K)=FO(2)*(B(I+1,J,K)+B(1-1,J,K)+B(1,J-1,K)+B(1,J+1,K)+ 
“ B(1,J,K+1)+B(1,J,K-1)-6. *B(1,J,K)) 
7 CONTINUE 


* B. THE TOP AND BOTTOM SURFACE NODES 


DO 8 I=2,26 
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DO 8 J=2,35 
BOUT(I,J,10)=FO(2)*(BCI+1,J,10)+B(I-1,J,10)+B(I,J+1,10)+ 
vs B(1I,J-1,10)+2*B(1I,J,9)+BI(2)*TINF-(BI(2)+6. )*B(I,J,10)) 
IF((I. GE. 9. AND. I. LE. 19. AND. J. GE. NC-1. AND. J. LE. NC+15) 
* | AND. . NOT. ((J. EQ. NC-1. OR. J. EQ. NC+15). AND. (I. EQ. 9. OR. I. EQ. 
* 19))) GOTO 8 
BOUT(I,J,1)=FO(2)*(B(I+1,J,1)+B(1I-1,J,1)+B(1,J+1,1)+B(I,J-1,1) 
* +2*B(I,J,2)+BI(2)*TINF-(B1I(2)+6. )*B(I,J,1)) 
8 CONTINUE 


faeces LHE EXTERIOR END FACES (COARSE INTERFACE) 


DO 9 I=2,26 
LTA=INT((I-. 5)/3)+7 
DO 9 K=2,9 
BOUT(I,1,K)=FO(2)*( BC I+1,1,K)+B(I-1,1,K)+B(1I,2,K)+B(1,1,K+1)+ 
x B(I,1,K-1)+. 5*A(IA,NB-1)-5. 5*B(I,1,K)) 
BOUT(I ,36,K)=FO(2)*(BCI+1, 36 ,K)+B(I-1,36,K)+B(I,35,K)+ 
B(1,36,K+1)+B(1,36,K-1)+. 5*A(IA,NB+12)-5. 5*B(1I,36,K)) 
9 CONTINUE 


* D. THE EXTERIOR SIDE FACES (COARSE INTERFACE) 


DO 10 J=2,35 
JA=INT((J-. 5)/3)+NB 
DO 10 K=2,9 
BOUT(1,J,K)=FO(2)*(B(2,J,K)+B(1,J+1,K)+B(1,J-1,K)+B(1,J,K+1)+ 
t B(1,J,K-1)+. 5%A(6,JA)-5. 5*B(1,J,K)) 
BOUT( 27,J,K)=FO(2)*(B(26 ,J,K)+B(27,J+1,K)+B(27,J-1,K)+ 
x B(27,J,K+1)+B(27,J,K-1)+. 5%A(16,JA)-5. 5*B(27,J,K)) 
10 CONTINUE 


eee 1HE EXTERIOR END EDGES (COARSE INTERFACE) 


DO 11 I=2,26 
IA=INT((1-.5)/3)+7 
BOUT(1,1,1)=FO( 2)*(B(I+1,1,1)+B(I-1,1,1)+B(I,2,1)+B(1,1,2)+ 
ve . SATA ,NB-1)+BI(2)*TINF-(4, 5+B1(2))*B(I,1,1)) 
BOUT(I,1, 10)=FO(2)* CECT TONPR Cie 1 on 4b 10)48Ct 9 )+ 
5#A(IA,NB-1)+BI(2)*TINF-(4. 5+BI(2) )*B(I,1,10)) 
BOUT(I ,36,1)=FO(2)**(BCI+1,36,1)+B(1-1,36,1)4+B(1,35,1)+B(1,36,2)+ 
ve ’5*A(IA,NB+12)+BI(2)*TINF-(4. 5+BI(2))*B(1,36,1)) 
BOUT(I ,36,10)=FO(2)**(B(I+1,36,10)+B(I-1,36,10)+B(1I,35,10)+ 
t B(1,36,9)+. 5*A( IA,NB+12)+BI(2)*TINF-(4. 5+BI(2))*B(I,36,10)) 
11 CONTINUE 


* FF. THE COARSE SIDE EDGES (COARSE INTERFACE) 
DOei Ze s=2), 32 


JA=INT((J-. 5)/3)+NB 
BOUT(1,J,1)=FO(2)*(B(2,J,1)+B(1,J+1,1)+B(1,J-1,1)+B(1,J,2)+ 


% .5%*A(6,JA)+BI(2)*TINF-(4. 5+B1(2))*B(1,J,1)) 
POUMm@oed 2 )=FO02)-(8( 26,0, 1)4B( 27 , J+), 1)48(27,0-1,1)4+8027,J,2)+ 
ve . 5*A( 16, JA)+BI(2)*TINF-(4. 5+BI(2))*B(27, J;1)) 


BOUT(1,J,10)=FO(2)**(B(2,J,10)+B(1,J+1,10)4+B(1,J-1,10)+B(1,J,9)+ 
% 5%A(6,JA)+BI(2)*TINF-(4. 5+BI(2))*B(1,J,10)) 
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BOUT(27,J,10)=FO(2)*(B( 26,J,10)+B( 27 , Jig 10) 4BC 27a = one 
* B(27,J,9)+. 5*AC16,JA)+BIC2)*TINF-(4. 5+B1(2)) B27 oe 
12 CONTINUE 


*% G. THE EXTERIOR CORNER EDGES (COARSE INTERFACE) 


DOM Suh —249 
BOUT(1,1,K)=FO(2)*(B(2,1,K)+B(1,2,K)+B(1,1,K+1)+B(1,1,K-1)+ 
* . 5%*(A(6,NB)+A(7,NB-1))-5. *B(1,1,K))- 
BOUT( 27, 1,K)=FO(2)*(B(26,1,K)+B(27,2,K)+B(27,1,K+1)+B(27,1,K-1)+ 
* -5%*(A(16,NB)+A(15,NB-1))-5. *B(27,1,K)) 
BOUT( 1,36,K)=FO(2)*(B(2,36,K)+B(1,35,K)+B( 1,36, K+1)+B(1,36,K-1)+ 
x .5%*(A(6,NB+11)+A(7,NB+12))-5. *B(1,36,K)) 
BOUT( 27, 36,K)=FO(2)*(B( 26, 36,K)+B( 27, 35,K)+B(27,36,K+1)+ 
Xe B(27,36,K-1)+. 5*(A( 16,NB+11)+A(15,NB+12) )-5. *B(27,36,K)) 
13 CONTINUE 


* H. THE EXTERIOR CORNERS (COARSE INTERFACE) 


BOUT(1,1,1)=FO(2)*(B(2,1,1)+B(1,2,1)+B(1,1,2)+. 5*(A(6,NB)+ 
ve A(7,NB-1))+BI(2)*TINF-(4. +B1I(2))*B(1,1,1)) 


Be 1, 10)=FOC2)*(BU2, I, 10 eae »10)+B(1, 1,9)+. 5*(AC6, NEOs 
A(7, NB- 1) )+B1(2)*TINF- (4. +BI(2))2 ‘BCL, 1, LO) 
“BOUT(1, Boye 1)=FO(2)* 1CBC2Z23605 1 )+B C135 7)eb ee 36 »2)+. 5*(AC6, NBs 


* A(7 ,NB+12))+BI(2)*TINF-(4. +BI(2))*B(1,36,1)) 
BOUDGI 3 cr 10)=FO(2)(B(2, 36,10)+B(1,35 10)+B(1, 36 59 )+. 5% 
Je (A(6, NB+11)+A(7, NB+12))+BI(2)* “TINE = (4; +B1(2))? ‘BC 1, 365 hore 


BOUT( 27,1, 1)=FO(2)*(B(26,1,1)+B(27,2,1)+B(27,1,2)+. 5*(A( 16 ,NB)+ 
* AC 15 ,NB-1))+BI(2)*TINF-(4. +BI(2))*B(27,1,1)) 
BOUT( 27,1, 10)=FO(2)**(B( 26,1, 10)+B(27,2, 10)+B(27,1 oo 
vs (A( 16, NB)+A(15 ,NB- 1))+BI(2)*TINF-(4. +BI(2))*B(27,1,10)) 
BOUT( 27,36, 1)=FO(2)%*(B( 26,36,1)+B(27,35,1)+B(27, 36 2) ieee 
ve (A(16, NB+11)+A(15, NB+12))+BI(2)*TINF-(4. +B1(2))**B(27,36,1)) 
BOUT( 27, 36, 10)=FO(2)**(B( 26, 36, 10)+B( 27,35, 10)+B(27,36,9)+. 5* 
t (AC 16 ,NB+11)+A(15 ,NB+12))+BI(2)*TINF-(4. +B1(2))*B(27,36,10)) 


* I. THE BOTTOM OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 20 I=10,18 
DO 20 J=NC,NC+14 
mee) 
DO 21 IC=3*(I-10)+1, 3*(1-10)+3 
DO 21 JC=3*(J-NC)+1, 3*(J-NC)+3 
21 TC=TC+C(IC,JC,8) 
BOUT(I ,J,4)=FO(2)*(B(I+1,J,4)+B(1-1,J3,4)+B(1,J+1,4)+B(1,J-1,4) 
“ +B(I,J,5)+TC/6. -6. 5*B(1,J,4)) 
20 CONTINUE 


* J. THE SIDE OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 22 J=NC,NC+14 
DO 22 K=2,3 
TC1=0. 
TC2=0. 
DO 23 JC=3*( J-NC)+1,3*(J-NC)+3 
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sc 


’ 
riy 


WOn23 KC=3"(K-2)+3,3*(K-2)45 
TTL, afc de 
23 He O24 C27 1G RC) 
BOUT(9,J,K)=FO(2)*(B(8,J,K)+B(9,J+1,K)+B(9,J-1,K)+B(9,J,K+1)+ 
* B(9,J,K-1)+TC1/6. -6. 5*B(9,J,K)) 
BOUT( 19, J,K)=FO(2)*(B(20,J,K)+B(19,J+1,K)+B(19,J-1,K)+ 
x B(19,J,K+1)+B(19,J,K-1)+TC2/6. -6. 5*B(19,J,K)) 
22 CONTINUE 


K. THE ENDS OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 24 I=10,18 
DO 24 K=2,3 
TC1=0. 
TC2=0. 
DO 25 IC=3*(I-10)+1,3*(I-10)+3 
BONS Gao ( K-27 ao 3% ( K-2) +5 
FOI=TCl+C( IG, 1, Ke) 
25 Reo —=TO2+G( IC, 45 KC) 
BOUT(I ,NC-1,K)=FO(2)**(B( I+1,NC-1,K)+B(I-1,NC-1,K)+B(I,NC-2,K) 
% +B(I,NC-1,K+1)+B(1I,NC-1,K-1)+TC1/6. -6. 5*B(I,NC-1,K)) 
BOUT( I, NC+15 ,K)=FO(2)*(B( I+1,NC+15,K)+B(1-1,NC+15 ,K)+B(I,NC+16,K 
“ )+B(1,NC+15 ,K+1)+B(I,NC+15 ,K-1)+TC2/6. -6. 5*B(I,NC+15 ,K)) 
24 CONTINUE 


i. THE SIDE EDGE OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 26 J=NC,NC+14 
TC1=0. 
TE2=0. 
DO 27 JC=3*( J-NC)+1,3*(J-NC)+3 
GoleTeire 1. Je. 1)+2*C01,JG,2) 
Ppmeeie2=10270( 27, JG, 1)+2*C( 27, JC, 2) 
BOUT(9 ,J,1)=FO(2)**(B(8,J,1)+B(9,J+1,1)+B(9,J-1,1)+B(9,J,2)+ 
ve TC1/6. +BI(2)*TINF-(5. 5+BI(2))*BC9,J,1)) 
BOUT(19,J, 1)=FO(2)**(B(20,J,1)+B( 19, J+1,1)+B(19,J-1,1)+B(19,J,2)+ 
x TC2/6. +B1(2)*TINF-(5. 5+BI(2))*B(19,J,1)) 
26 CONTINUE 


feelin END EDGE OF THE EXCLUSION ZONE (FINE INTERFACE) 


DO 28 I=10,18 
TC1=0. 
TC2=0. 
DO 29 IC=3*(I-10)+1,3*(1-10)+3 
Aelame ie te. 1. 1)+27C( iG. 1.2) 
Pon 9 fG2=1624C( 1C ,45,1)+2*CC1C,45,2) 
BOUT( I ,NC-1, 1)=FO(2)*(B(I+1,NC-1, 1)+B(I-1,NC-1, 1)+B(I,NC-2,1) 
x +B(I,NC-1,2)+TC1/6. +B1(2)*TINF-(5. 5+B1(2))*B(I,NC-1,1)) 
BOUT( I ,NC+15 , 1)=FO(2)**(B(14+1,NC+15,1)+B(1I-1,NC+15,1)+B(1I,NC+16, 1) 
x +B(I,NC+15,2)+TC2/6. +BI(2)*TINF-(5. 5+BI(2))*B(I,NC+15,1)) 
28 CONTINUE 


RETURN 
END 


14] 


4. Source Listing of Subroutine CORC 


SUBROUTINE CORC(A,B,AOUT) 
THIS SUBROUTINE HAS BEEN MODIFIED TO USE EXTENDED ARRAYS IN THE J 
DIRECTION TO ALLOW THE STUDY OF COOLING RATES. 


ao 


REAL A,B,AOUT 
DIMENSION A( 21,36) ,AOUT( 21,36) ,B(27,36,10),FO(3),BI(4) 
COMMON NB,NC,YARC,SUM,TINF,BI,FO 


* A. THE CORNERS 


AOUT(1,1)=0. 
AOUT(1,36)=0. 
AOUT( 21,1)=0. 
AOUT( 21, 36)=0. 


*) By ene eENDS 


DO 1 I=2,20 
AOUT(I,1)=FO(1)**(A( 1+1, 1)+A( 1-1, 1)+A(1,2)+TINF*BI(1)-(3+BI(1))* 
. BC Us), 
AOUT(I ,36)=FO(1)*(ACI+1, 36)+A(I-1, 36)+A(1, 35) +TINF*BI(1)-(3+BI(1) 
< Vee ele, oO) 
1 CONTINUE 


o!, 


*  C. THE SIDES 


DO 2 J=2,35 
AOUT(1,J)=FO(1)*(A(1,J+1)+A(1,J-1)+A(2,J)+TINF*BI(1)-(3+BI(1))* 
si Ol) ) 
AOUT( 21,J)=FO( 1)*(A(21,J+1)+A(21,J-1)+A( 20, J)+TINF*BI(1)-(3+ 
ve I Ca 
2 CONTINUE 


~ Do THE MiDSeOmNTs 


NMIN=NB-1 
NMAX=NB+12 
DO 3 I=2,20 
DO 3 J=2,35 
IF (J. GE. NMIN. AND. J. LE. NMAX. AND. I. GE. 7. AND. I. LE. 15) GOTO 3 
IF (J. GE. NB. AND. J. LE. (NB+9). AND. (I. EQ. 6. OR. I. EQ. 16)) GOTO 3 
AOUT(I ,J)=FO( 1)*( AC 1+1,J)+A(1I-1,J)+A(1, J+1)+A(1,J-1)+BIC1)*TINF 
2 (G7 BIC 1) ACen 
3 CONTINUE 


% E. THE MEDIUM ENDS TRANSITION 


pe) ay 1S 

TA=0. 

TB=0. 

OMS 1B=(l-7 73416 (1-6) *3 
TASTES oly, oy 
Pb, 1.10) 2.474 
TB=B(1B,36,1)/2.+TB 
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TB=B(1B,36,10)/2. +TB 
DO 5 K=2,9 
TA=B(IB,1,K)+TA 

5 TB=B(1B,36,K)+TB 


AOUTC I ,NMIN)=FO( 1)*(ACI+1,NMIN)+AC I-1,NMIN)+AC1I,NMIN-1)+TA/18. 


* BI(1)*TINF-(4. 5+BI(1))*A(I ,NMIN) ) 


AOUT( I ,NMAX)=FO( 1)7*( A( I+1, NMAX)+A( I-1,NMAX)+A( I ,NMAX+1)+TB/18. 


* BI(1)*TINF-(4. 5+B1(1))**A(I , NMAX) ) 
4. CONTINUE 


* F. THE MEDIUM SIDE TRANSITION 


DO 6 J=NB,NB+11 

TA=0. 

TB=0. 

DO 7 JB=(J-NB)*3+1,(J-NMIN)*3 
TA=B(1,JB,1)/2.+TA 
TA=B(1,JB,10)/2.+TA 
TB=B(27,JB,1)/2.+TB 
TB=B(27,JB,10)/2.+TB 
DO 7 K=2,9 

TA=B(1,JB,K)+TA 
; TB=B( 27 ,JB,K)+TB 
AOUT( 6 ,J)=FO(1)*(A(5 ,J)+A(6,J-1)+A(6 , J41)+TA/18. + 


* BI(1)*TINF-(4. 5+BI(1))*AC6,J)) 
HPOUTGIO. JNO, 1) CAC? , J) +AC 16 ,J-1)+AC16,J+1)+1B/18.+ 
* BIC 1)*TINF-(4. 5+BI(1))*AC16,J)) 
6 CONTINUE 
RETURN 
END 


5. Source Listing of Subroutine RATES 


SenROUDINE RATES(B,C,CK,CT,CT2,NB,NC,TIME) 


a 


a 


* THIS PROGRAM DETERMINES THE COOLING RATES WHILE RUNNING THE PROGRAM 


WELDC. 
DIMENSION C(27,45,8),CR(200,14,3,4) ,CT(200,14,3),CT2( 200,14, 3) 
* B(27,36,10) 
mMOGIOAL, CT,CT2 
LOCAL=9**NB+3**NC-56 
C CHECK THE TEMPERATURES IN THE FINE ZONE 
DO 1 I=1,14 
J=0 
DO 1 JP=LOCAL, LOCAL+44 
J=J+1 
C CHECK TO SEE IF ABOVE THE UPPER LIMIT 
Mreceeted. 1). GE, (825.0)) THEN 
CT(JP,1,1)=. TRUE. 
ENDIF 
mam@eciedi.1).CE.(675.0)) THEN 
CT(JP,1,2)=. TRUE. 
ENDIF 
if 066@0 9 11)-GE.(1075.0).) THEN 
CT(JP,1,3)=. TRUE. 
ENDIF 


C CHECK TO SEE IF PREVIOUSLY ABOVE UPPER LIMIT BUT NOW BELOW. IF S80 


C SAVE TIME AND TEMPERATURE. 
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IF (C(I,J,1). LT. (825. ). AND. CT(JP,1,1). AND. CoNOT, Cl2( JP sien 
* THEN 
CR(JP,1I,1,1)=C(I,J,1) 
CR(JP,1I,1,2)=TIME 
CT2(JP,1,1)= TRUE 
ENDIF 
IF (C(I,J,1). LT. (673. 0). AND. CT( JP,1,2). AND. (, NODWCT2( JPalenm 
* THEN 
CROIE 1 22) Olean 
CR(JP,1,2,2)=TIME 
CT2(JP,1,2)=. TRUE. 
ENDIF 
IF (C(I,J,1). LT. (1073. 0). AND. CT( JP,1,3). AND. (. NOT) CT2( Jpmieaem 
* THEN 
CRelieelen sen) Clee gi aT) 
CR(JP,1,3,2)=TIME 
CT2(JP,1,3)=. TRUE. 
ENDIF 
C IF TEMPERATURE IS LESS THAN THE LOWER LIMIT, SAVE THE TEMPERATURE AND 
C THE TIME. 
IF (CT2(JP,1,1). AND. (C(1,J,1). LT. 795. ). AND. CT(JP,I,1)) THEN 
CT(JP,1I,1)=. FALSE. 
GR( JE, Tien =o 1 ed) 
GR(JP,1,1,4)=TIME 
ENDIF 
IF (CT2(JP,1,2). AND. CCC1.J.1). LT, 423) AND. GT( JP. 1,2) ene 
CT(JP,1,2)=. FALSE. 
CRIP 1.2. 3)=Coneie a) 
CR(JP,1,2,4)=TIME 
ENDIF 
IF (CT2(JP,I,3). AND. (C(I,J,1). LT. 773. ). AND. CT(JP,1,3)) THEN 
CT( JP,1,3)=. FALSE. 
CR(JP,1I,3,3)=C(I,J,1) 
CR(JP,1,3,4)=TIME 
ENDIF 
1 CONTINUE 
C CHECK THE TEMPERATURES IN THE MEDIUM ZONE 
DO 2 I=10,14 
IP=(1-9)*3-1 
DO 2 J=1,NC-1 
JP=9*NB+J*3-55 
ih (eee aay 2 
C CHECK TO SEE IF ABOVE THE UPPER LIMIT 
IF (B(I,J,1).GE.(825.0)) THEN 
CT(JP,IP,1)=. TRUE. 
ENDIF 
IF (B(I,J,1).GE.(675.0)) THEN 
Cie Beg) — ee 
ENDIF 
IF (B(1I,J,1).GE.(1075.0)) THEN 
CT(JP,IP,3)=. TRUE. 
ENDIF 
C CHECK TO SEE IF PREVIOUSLY ABOVE UPPER LIMIT BUT NOW BELOW. IF SO 
C SAVE TIME AND TEMPERATURE. 
IF (B(I,J,1). LT. (825. ). AND. CT(JP,IP,1). AND. (. NOT. CT2(JP,IP,1))) 
*© THEN 
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eRqe. (Piel j=BC 1 1) 
CR(JP,IP,1,2)=TIME 
CT2(JP,IP,1)=. TRUE. 
ENDIF 
feet e( 1.) oninC67 sO mAnbeGiG Ie. IPe2y AND. ¢. NOT.Gm2@JP,1P,2))) 
* THEN 
CR(JP,IP,2,1)=B(1,J,1) 
CR(JP,1P,2,2)=TIME 
G02 (JP. 1P.2)=. TRUE. 
ENDIF 
IF (B(1,J,1). LT. (1073. ). AND. CT( JP, IP,3). AND. (. NOT. CT2(JP,IP,3))) 
* THEN 
CR(JP,IP,3,1)=B(1,J,1) 
CR(JP,IP,3,2)=TIME 
CT2(JP,IP,3)=. TRUE. 
ENDIF 
C IF TEMPERATURE IS LESS THAN THE LOWER LIMIT, SAVE THE TEMPERATURE AND 
C THE TIME. 
IF (CT2(JP,IP,1). AND. (B(1,J,1).LT. 795. ). AND. CT(JP,IP,1)) THEN 
CT(JP,IP,1)=. FALSE. 
CR(JP,IP,1,3)=B(1,J,1) 
CR(JP,IP,1,4)=TIME 
ENDIF 
iaeecon2( JPe lee AND. (8(1.J,1). LT. 623. ). AND. CICJP, 1P.2))) THEN 
CT( JP, IP,2)=. FALSE. 
GRGr 1P.2.3)-n( 1) 
CR(JP,IP,2,4)=TIME 
ENDIF 
IF (CT2(JP,IP,3). AND. (B(1,J,1). LT. 773. ). AND. CT(JP,IP,3)) THEN 
CT( JP, IP,3)=. FALSE. 
GROIP IP, 3. Sy—nel. J, 1) 
CR(JP,IP,3,4)=TIME 
ENDIF 
2 CONTINUE 
RETURN 
END 


6. Source Listing of Program RATEOUT 


PROGRAM ROUT 
C USED TO OUTPUT COOLING RATES CALCULATED BY THE PROGRAM PLOTC 
fee -7-SS8 BY ROBERT ULE 


MOGTGARTONC200, 14,3) ,Gh2( 2008 1453) 
REAL*4 CR(200,14,3,4) ,CRATE( 200, 14,3),RAD(200, 14,3) ,Y(200,14,3) 
DATA GCRATE/8400*0. /,RAD/8400*0. /,Y/8400*0. / 


OPEN (1,FILE='RATE' ,STATUS='OLD' , FORM='UNFORMATTED’ ) 
OPEN (2,FILE='ROUT' ,STATUS='NEW' ) 

READ( eGR Cr ,CT2 

CLOSE (1) 


C CALCULATE COOLING RATE IF COOLED THROUGH THE TRANSITION ZONE 
fepobloE, INDICATE WHERE IT TEMPERATURE CYCLE THAT POINT Is. 


DO 10 J=1,200 
OG GO) eth ibe 
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10 


DO slew 

IF ((.NOT. CT(J,1,K)). AND. (. NOT. CT2(J,1,K))) CRATE(J,1,K)=-3. 
IF (CT(J,1,K). AND. (. NOT. CT2(J,1,K))) CRATE(J, 1,k)=-2, 
IF (CT(J,1,4) cAND.CI2C) (1 .K))) CRATE GI) eae 

IF ((.NOT. CT(J,1,K)). AND. CT2(J,1,K)) THEN 
CRATE(J,1,K)=(CR(J,1,K,1)-CR(J,1,K,3))/(CR(J,1,K,4) 

< BO WS 5 23))) 

YARC=34, +22 “(CROJ Kk, 2)CR0 91K ae 
Y(J,1,K)=YARC-J 
RAD(J,1,K)=SQRT(Y(J,1,K)**2+(14-1)**2) 

ENDIF 

CONTINUE 


C OUTPUT THE COOLING RATE DATA FOR THE THREE DIFFERENT CRITERIA AND 
C AT EACH POINT 


20 


30 


40 
60 


70 


80 


GOTO 60 
WRITE(2, (1X; 4@4,,/)) _) 
*' COOLING RATE OUTPUT KEY: ', 
*'-3 = STILL BELOW TRANSITION TEMPERATURE UPPER LIMIT’, 
*'-2 = ABOVE TRANSITION TEMPERATURE UPPER LIMIT', 
*'-1 = COOLING BETWEEN UPPER AND LOWER TRANSITION LIMITS' 
DO 20eJ=207 A100 

WRITE(2,'(/A11,13,/)')'Y POSITION ',J 

WRIDE( 2, CA yee POS Y R MORRIS' 

Om Onl — 9 aly 

WRITE(2, (iXalom2k 266s 2.2 Gneeye le 


* VC Jol eR Ames el 1 CRATE Crit) 
CONTINUE 
DO 30 J—30 eo 
WRITDOC2. (/All, f3./)) lox eeOSMMION J 
WRITE(2,'(A)') “X EOS Y R MATERIAL’ 
DO 30 I=8,14 
WRITE( 23 7 Gixiige 2X, OF 6.2. bil 2a.) ole 
ve VCJ, la) RADCd. 2b 2 CRIP Gens. 
CONTINUE 
DO 40 J=30,100 
WRITE( 2, (/Al1.13,/) )'Y ROSMaG head 
WRITE(2,'(A)') 'X POS Y R HYDROGEN’ 
DOn40 1=8 14 
WEMME G2pemGlx 1 3.02%. 21 O22 sh eee eee 
* M¢Jyias) ,RADCJ,1,3)7 CRAEECS, LoD 
CONTINUE 
OPEN(10,FILE='PLOT1' ,STATUS='NEW' ) 


DO 70 J=30,96 

IF (CRATE(J,14,1).GT.0) WRITE(10,'(2F12.4)') FLOAT(J-36), 
* (Rada), 14,1) 

CONTINUE 

OPEN(11,FILE='PLOT2' , STATUS=' NEW’ ) 

DO 80 J=30,96 

IF (CRATE(J,13,1).GT.0) WRITE(11, (2F12.4). ep eoAl(s=4e 
eICRADE CI) «borat ) 

CONTINUE 

OPEN(12,FILE='PLOT3' , STATUS=' NEW’ ) 

DO 90 J=30,96 

IF (CRATE(J,12,1).GT.0) WRITE( 12) (2FlZ3 2) Fle a@-369e 
* ORATEC J, 2 19 
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90 CONTINUE 
OPEN(13,FILE=' PLOT4' , STATUS='NEW' ) 
DO 95 J=30,96 
Pe eCGRATECS. 11.1) CLO aWrimie@ls. (2Pi2.4)°) FLOAT(J-36) , 
CGRAnEC I. lie 1) 
95 CONTINUE 
OPEN(14,FILE='PLOTS' , STATUS=' NEW’ ) 
DO 96 J=30,96 
BEMCORATEC J, 10,1).G0.0) WRITECla, (2h12. 4) ) FROAT(J-36), 
Meer ATE( J ,10,1) 
96 CONTINUE 
OPEN(15,FILE=' PLOT6' ,STATUS=' NEW’ ) 
DO 97 J=30,96 
IF (CRATE(J,9,1).GT.0) WRITE(15,'(2F12.4)') FLOAT(J-36), 
* CRATE(J,9,1) 
97 CONTINUE 
END 
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